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FREQUENCY  DOMAIN  INTERPRETATION  OF 
ITERATIVE  DECONVOLUTION  TECHNIQUES 


INTRODUCTION : 

Estimating  the  input  to  a  linear  systems  when  the  output  (corrupted 
by  noise)  and  the  system  function  are  known,  is  called  the  deconvolution 
problem.  Many  approaches  to  the  solution  of  the  problem  involve  iterative 
time  domain  techniques.  Here,  we  provide  a  frequency  domain  interpretation 
of  these  techniques. 

Formulation  of  the  Problem 

The  problem  may  be  viewed  mathematically  as  follows: 

Define: 


Input 

Kt) 

Output 

0(t) 

Instrument  function 

h(t) 

Measured  noise 

n(t) 

Since  the  instrument  is  a  linear  system,  I(t)  and  0(t)  are  related  by  the 
convolution  integral  (*  denotes  convolution) . 

O(t)  «  I(t)  *  h(t)  +  n(t) 


or 

0(t)  KO  h(t-x)  dT  +  n(t) 

In  the  problem  we  are  addressing,  0(t)  and  h(t)  are  known  (measured  and 
we  wish  to  estimate  I(t). 

Solution  Methods 

The  reference  papers  provide  a  summary  of  some  of  the  approaches  that  have 
been  employed  to  solve  this  problem. 
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Most  of  those  techniques  employ  iterative  methods. 

Historically  the  noise  term  has  been  (formally)  ignored  and  the 
deconvolution  has  been  viewed  as  a  Fredholm  integral  equation  of 
the  first  kind.  The  iterative  techniques  which  have  been  employed 
use  the  concept  of  successive  substitutions.  (1,2) 


The  basic  types  of  algorithm  in  these  iterative  methods  are  a 
variation  of  the  following:  with  0(t)  the  measured  output, 

I1  (t)  =  0(t)  +  [  0(t)  -  0(t)  *  h(t)  ] 

I  (t)  -  I  .(t)  +  [  0(t)  -  I  ,(t)  *  h(t)  ], 
n  n-i  n-o. 

f inally 

A 

I(t)  =  1  (t)  for  some  predetermined  n 
n 

A 

(  I,  denotes  the  estimated  input) 
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Minimum  Mean  Squared  Error  Solution  (MMSE) 

Viewing  this  problem  as  similar  to  the  system  identification  problem 
where  h(t)  is  unknown  [12,  13],  consider  deriving  the  MMSE  estimate 
for  I(t).  The  MMSE  estimate  would  have  the  interpretation  of  having 
the  smallest  variance  of  any  possible  estimate  of  I(t). 


Thus  from  Eq.  (1)  we  wish  to  choose  I(t)  to  minimize  the  expected 

2 

value  (statistically)  of  £  where: 


£  =  E  |N(t)}  =  E  {  0(t)  -  J  I(T)  h(t-T)  d  T  | 

That  is  find  I(t)  such  that 

Min  £2  =  Min  E  {  N(t)2} 

(2)  =  Min  E  {  \  0(t)  -f  I(t)  h(t-x)  dT  }  2  } 

where 

E  denotes  the  mathematical  expectation  (statistical  average) .  Under  the 
assumption  that  N(t)  is  not  correlated  with  h(t),  it  is  well  known  (as  a 
consequence  of  the  orthogonality  principle)  that  the  solution  to  F.q.  (2) 
leads  to  the  Wiener-Hopf  equation. 


(3) 

where 


©oh(t)  =  f  I(T  )  9h  (t-Tj  dT 


0,  is  the  autocorrelation  of  the  system  function  h(t) 
h 

and  &  ,  is  the  cross  correlation  of  the  output  and  system  functions, 
oh 

0(t)  and  h(t)  respectively. 


This  equation  may  be  solved  for  I(t)  using  Fourier  transform  techniques. 
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Letting 

F  0  oh  <*> 

f  e  h2  (t) 


®oh(b,) 
®h2  ((0) 


and 

F  jl(t) 

we  obtain 


If  (o>)  , 


(4) 

IF  (to)  = 


(assuming  0,2 
h 


^oh  ^ 

®h2  (co) 

(co)  is  not  zero  over  the  range  of  to  considered) 


Finally  taking  the  inverse  Fourier  transform 


(5)  * 

I  (t) 


^oh 

v  <«> 


] 


Note  that  the  FFT  may  be  employed  to  perform  the  required  calculations. 

If  there  is  no  noise  i.e.  N(t)  =  0,  (5)  will  yield  the  exact  solutions 

to  (1). 

When  noise  is  present  some  smoothing  of  the  spectral  estimate  ,  (to)  will 

oh 

be  required  to  produce  a  stable  estimate. 
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There  is  a  basic  trade-off  between  increased  resolution  of  peaks 
in  I(t)  and  noise  performance.  That  is,  in  the  presence  of  noise, 
erroneous  noise  spikes  appear  in  I(t)  as  increased  resolution  is 
demanded.  In  fact,  with  no  spectral  smoothing  [highest  resolution 
possible]  virtually  no  noise  can  be  tolerated.  This  is  due  basically 
to  the  fact  that,  without  smoothing,  the  MMSE  criterion  minimizes  the 
square  error  integrated  over  the  signal  space.  Hence,  MMSE  is  some¬ 
what  insensitive  to  noise  spikes. 

Iterative  Method  as  a  Special  Case  of  the  MMSE  Approach 

Following  the  presentations  of  [2]  and  [8],  one  can  show  that  the 
iterative  solution  (after  n  iterations)  can  be  expressed  as  follows: 

In(t)  =  (n+DOCt)-^1)  h(t)  *  0(t) 

+  (H*1)  h(t)  *  h(t)  *  0(t) 

+ - +  h(t)  * - .  *  h(t)  *  0(t) 

n  times 


Taking  the  Fourier  transform  of  both  sides  and  collecting  terms  we  have 

(n+1) 


IF(w) 


n 

+  E- 
m=l 

(-1)" 

(n+1) 

'm+1' 

m 

hf 

(u>) 

X) 

=  E 

TT1“0 

(-i)“ 

(n+1) 

'm+1' 

H  “ 
t 

(CO) 

* 

=  hf 

(co)  .  oF 

(w) 

1 

1  H: 

F(to)  1 2 

• 

nr 

(-D"  <2J>  H^1  (co) 


Unsmoothed  MMSE 


Smoothing  Function 
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Using 


®oh  (U) 


H*p  (to)  ,  O Tiu>) 

J  Hy  (to)|  2  "  1  Oh2  (OJ)  J 

in  combination  with  (4),  (5)  and  this,  last  equation  shows  that  the 
iterative  method  may  be  viewed  as  the  MMSE  method  with  the  spectral 


smoothing 

function 

given  by: 

S  (to) 

□ 

=  X 

m=o 

(-D" 

n+l 

hn+l' 

m+l  ,  . 

(to) 

or 

nil 

Jl 

(-1)"'1 

(n^) 

■s 

(co) 

This  is  depicted  in  Figure  1  below. 

Parallelling  the  derivations  (12]  and  [13]  one  can  establish  (assuming 
Gaussian  noise  N(t))  that 

A 

Ip  (to)  is  an  F  distributed  random  variable  where  the 
degrees  of  freedom  dpend  on  the  smoothing  function  S  (to)  . 


0(t) 

Ideal 

Smoothing 

\ 

Inverse 

Filter 

S 

Filter 

7 

S(to) 

^Kt) 


To  improve  noise  performance) 


Frequency  Domain  Interpretation  of  Iterative  Deconvolution  Techniques. 


FIGURE  1 
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AN  ESTIMATE  OF  THE  EFFECTS  OF  TEMPORAL 
_  BACKGROUND  ON  SYSTEM  PERFORMANCE 


1.  INTRODUCTION 

In  order  to  evaluate  the  performance  of  an  infrared  sensor,  the  system 
designer  needs  to  knoy  the  spatial  and  temporal  characteristics  of  the 
background.  With  these  data  and  a  knowledge  of  platform  instabilities  he 
can  systematically  design  a  sensor  and  predict  its  sensitivity,  false  alarm 
rate  and  detection  probability.  When  one  considers  a  mosaic  staring  sensor, 
spatial  and  temporal  variations  in  the  background  affect  the  performance  in 
two  ways.  First,  platform  instabilities  can  cause  radiance  variations  from 
the  spatially  structured  terrain  and  second  time  variant  changes  in  the 
background  will  cause  a  response  in  the  sensor. 

Over  the  past  several  years,  techniques  for  the  performance  analysis 
of  advanced  surveillance  systems  have  been  developed.  Much  of  the  work  has 
been  limited  however  to  analysis  of  spatial  variations  in  the  background 
obtained  from  various  measurement  programs.  Among  these  are: 

.  RM-19  Satellite  Experiment  (Lockheed)  -  The  RM-19  experiment  was 

designed  to  obtain  background  measurements  from  a  low  altitude, 
high-inclination,  circular  orbit,  in  a  number  of  spectral  bands 
near  2,  A  and  6  ym.  Spectral,  spatial  and  exceedance  data  were 
obtained  and  represent  one  of  the  major  sources  of  background  data 
from  space. 

.  BMP  -  CAMP  Experiment  (Lockheed  U-2)  -  The  Lockheed  U-2  program  was 
designed  to  measure  spatial  characteristics  of  various  bakground 
sources  in  the  blue  and  red  spike  region  of  Am  as  well  as  at 
various  wavelengths  in  the  LWIR  region  14  ym. 
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Valuable  PSD  data  have  been  made  available  through  this  program,  which 
shows  the  dependence  on  terrain  and  cloud  types  as  well  as  the  varia¬ 
tion  of  PSD  from  MWIR  to  LWIR. 

.  KC-135  flights  (AFGL)  -  High  resolution  spectra:  to  =  4  cm  ^  and  ba  = 

1  cm  ^  of  several  classes  of  terrain  and  cloud  backgrounds  have  been 
measured  from  KC-135  aircraft.  These  data  have  been  used  to  establish 
bench  marks  to  calibrate  the  emission  and  reflection  models  from 
terrain  and  clouds  as  a  function  of  solar  zenith  and  azimuth  angles. 

.  U-2  flights  (A.D.  Little)  -  Additional  cloud  spectral  data  are  avail¬ 

able  from  the  A.D.  Little  U-2  flights.  These  consist  of  a  series  of 
observations  of  clouds  of  varying  types,  and  over  a  range  of  viewing 
conditions  in  the  2  ym  region.  These  data  are  useful  to  provide  a 
range  of  variation  to  be  expected  in  cloud  spectra. 

.  ERIM  flights  -  High  spatial  resolution  data  was  obtained  over  desert 
and  mountainous  terrain,  over  rural  and  urban  terrain  and  over  water. 
Spatial  resolution  was  two  to  three  feet.  Data  was  collected  in  six 
selected  bands  from  2  to  5.7  ym  as  well  as  the  spectral  band  from  2 
to  11.4  ym.  The  primary  interest  in  these  data  is  its  two  dimensional 
nature.  Unfortunately,  two  dimensional  data  are  not  presented,  the 
presentation  being  in-track  and  cross-track  power  spectra. 

All  of  these  programs  were  concerned  with  spectral  and  spatial  varia¬ 
tions  of  terrain  and  cloud  background?  Although  there  has  been  work  on 
modelling  temporal  backgrounds,^  temporal  background  data  has  not  been 
available.  It  is  expected  that  a  strong  source  of  temporal  fluctuations 
would  be  moving  clouds.  In  the  absence  of  experimental  data,  a  random 
process  model  was  postulated  to  obtain  power  spectral  density  estimates. 

It  is  the  purpose  of  this  paper  to  present  this  model  and  to  estimate  the 
effect  of  this  clutter  on  system  performance. 
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2.  LWIR  CLOUD  TEMPORAL  CLUTTER  MODEL 


Consider  a  broken  cloud  distribution  moving  past  a  starting  sensor. 

If  we  assume  the  cloud  is  opaque,  the  radiance  observed  at  the  sensor 

represents  the  radiation  from  a  blackbody  at  the  temperature  of  the 

ambient  air  at  the  cloud  altitude.  When  the  atmosphere  is  clear,  the 

sensor  would  observe  the  warm  terrain.  To  illustrate  this  situation 

a  hypothetical  trace  is  shown  in  Figure  1.  This  illustrates  the  radiance 

fluctuations  that  a  downward  looking  sensor  would  observe  generally  at  a 

wavelength  of  9  ym  for  clouds  at  3  or  6  km.  It  is  reasonable  to  assume 

2 

that  the  random  arrival  of  clouds  may  be  modelled  as  a  Poisson  process 
with  arrival  parameter  adjusted  for  an  assumed  percentage  of  cloud  cover. 
For  the  calculations  that  follow  a  50%  cloud  cover  is  assumed.  A  typical 
sample  of  Poisson  arrivals  is  shown  in  Figure  2.  The  effect  of  a  cloud  on' 
the  measured  radiance  depends  on  the  wavelength  and  the  altitude,  speed, 
size  and  shape  of  a  cloud  passing  through  the  field  of  view.  For  the 
purpose  of  this  paper,  we  have  assumed  a  sensor  operating  at  9  viewing 
terrain  with  broken  clouds  at  6  km  altitude.  The  idea  is  to  allow  clouds 
to  enter  the  f ield-of-view,  reach  a  maximum  image  of  radiance  and  then 
exit.  The  effect  of  a  cloud  on  the  radiance  trace  was  assumed  to  take 
one  of  the  forms  shown  in  figure  3.  Obviously,  many  shapes  may  be  assumed 
but  it  is  felt  that  those  shown  would  be  representative  of  the  problem 
allowing  for  an  upper  limit  of  fast  entry  and  exit  as  well  as  for  a  smooth 
entry  and  exit,  thus  bounding  the  probable  spectral  roll-off  that  could 
occur  with  an  actual  cloud. 

The  time  required  for  entry  and  exit  and  the  duration  of  time  the 
cloud  is  in  the  f ield-of-view  depends  on  the  size  and  speed  of  a  cloud. 

The  radiance  waveform  (figure  3)  may  be  coupled  with  the  Poisson  arrival 
process  by  convolving  the  shape  with  the  Poisson  impulses.  This  yields 
what  is  commonly  called  a  "shot  noicc"  process  model. 
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If  we  denote  the  cloud  shape  trace  (Figure  4)  by  h(t) ,  the  assumed 
process  may  be  represented  as2 


h(t)  =  E  h(t  -  t.) 
i— «o 


(1) 


where  t^  are  Poisson  random  time  arrivals.  The  statistical  autocovariance 
function  for  the  process  g(t)  is  obtained  by  convolving  the  process  with 


itself.  Thus: 


/ 


C(T)  =  \  /  h(a)h(T  +  a  )  da 


(2) 


The  power  spectrum,  the  fourier  transform  of  C(T)  can  be  obtained  from  the 
impulse  response  h(t)  and  is  given  by 

S(co>  -  \  |  H(co>f  2  (3) 

The  fourier  transform  of  h(t)  is  given  as 

H(co)  =  F[h(t)  ]  (4) 

For  the  waveforms  of  Figure  3,  the  corresponding  H(co)  are  as  follows: 

-  (5)  reference  3 


(a) 

H  (co) 

=  ATo 

(b) 

H  (co) 

2A 

b-a 

(c) 

H  (co) 

=  A 

sin  co  T^/2 


"  V2 

cos  a  co  -  cos  b  co 


(6)  reference  3 


<0 


(1  -  2a) 


T  T  /2 

0  sin  co  O' 


Tj  2 

03  0 


sin  a>  V2  (l-2a) 
co  Tq/2  (l-2a) 


-  cos  a  co  TQ/2 


.  sin  co  TQ/2  (1  -  a) 


03 


2 

03  ~ 


4tx 


(7)  reference  4 

2 


2  2 
4a2  Tq2 
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3. 


CALCULATED  RESULTS 


For  the  results  that  follow,  footprints  of  200  m  and  1  km  were 

chosen  as  a  representative.  For  the  200  m  footprint  we  arbitrarily 

assumed  cloud  sizes  of  200  m,  500  n,  1  ka  and  5  km.  For  the  1  km 

footprint  we  assumed  cloud  sizes  of  1  km,  5  km  and  10  km.  We  assumed 

cloud  speeds  of  10  and  30  m/sec  which  are  not  unreasonable  speeds 

based  on  data  shown  in  Figure  4."*  As  indicated  earlier,  a  50%  cloud 

cover  was  assumed.  Cloud  size  and  velocity  were  separately  evaluated 

so  that  14  separate  cases  were  treated  for  both  the  simple  trapezoidal 

waveform  Equation  (6)  and  cosine  trapezoid  Equation  (7).  The  power 

spectrum  for  each  case  was  calculated.  Selected  PSD  are  presented  in 

the  following  figures  to  show  the  effects  of  cloud  size,  speed  and  type 

of  waveform.  All  clouds  are  assumed  to  be  at  6  km  altitude.  For 

clouds  at  3  km  the  results  should  be  divided  by  2.  Figure  5  shows  the 

PSD*  of  the  random  arrival  of  200  m  clouds  moving  at  a  speed  of  10  m/sec 

assuming  a  cosine  trapezoidal  waveform.  The  first  null  occurs  at 

approximately  0.05  Hz  with  a  roll  off  of  approximately  f  The  same 

situation  is  shown  in  Figure  6  except  for  a  trapezoidal  waveform.  The 

first  null  occurs  at  approximately  the  same  frequency.  However,  the 

-3.5 

roll  off  is  now  changed  to  f  .  The  effect  of  the  increase  in  velocity 

is  shown  in  Figure  7  for  the  same  conditions  shown  in  Figure  5.  The 

velocity  of  the  clouds  are  now  30  m/sec.  The  first  null  has  shifted  to 

-4 

a  higher  frequency  approximately  o.l5  Hz  and  the  roll  off  is  now  f 

The  situation  for  larger  footprint  is  shown  in  the  next  two  figures. 

Figure  8  shows  the  PSD  of  1  km  clouds  moving  at  speeds  of  10  m/sec  across 

a  1  km  footprint.  The  first  null  now  occurs  at  0.01  Hz  with  a  roll  off 

of  approximately  f  Figure  9  shows  the  PSD  of  10  km  clouds  moving  at 

speeds  of  30  m/ sec  across  a  1  km  footprint.  The  first  null  occurs  at  a 
frequency  of  approximately  0.002  Hz  with  a  roll  off  of  f 

^Valley,  S.L.,  Editor,  Handbook  of  Geophysics  and  Space  Environments ,1965 
*The  PSDs  should  read  log  PSD  in  the  figures. 
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The  rectangular  waveform  of  Figure  3  (a)  is  not  considered  realistic 

since  clouds  will  not  enter  or  exit  that  abruptly.  However,  this  wave- 

-2 

form  serves  as  an  upper  bound  on  the  high  frequency  roll  off  rate  f  .  The 
first  null  of  the  spectra  depends  on  the  duration  of  the  pulses  y2/T.  This 
is  consistent  with  the  results  shown  in  Figures  5  through  9.  Thus,  the  speed 
and  size  of  the  cloud  dictate  the  bandwidth  of  the  spectra. 

The  RMS  clutter  leakage  was  obtained  for  each  case  by  filtering  through 

first,  sencond,  and  third  differencing  algorithms.  For  200  m  footprints,  a 

frame  time-of  1  sec  was  used.  For  km  footprint,  a  frame  time  of  5  sec  was 

used.  Results  for  third  differencing  are  plotted  in  Figure  10.  In  the 

figure,  the  clutter  leakage  is  shown  vs.  cloud  size  with  cloud  speed  and 

footprint  as  parameters.  For  all  the  cases  considered,  the  clutter  leakage 

2 

varies  over  the  range  from  10-100  kw/cm  sr  ym. 

For  a  more  realistic  situation  one  can  combine  a  number  of  different 
cloud  sizes  moving  at  a  constant  speed  at  a  given  altitude  or  combine  clouds 
at  different  altitudes  with  different  speeds.  Assuming  they  are  independent 
the  spectra  would  add  since  the  cross  covariances  are  zero.  This  results 
in  a  spectra  model  given  by 

N  2 

S(uj)  =  2  PA.  H(o>)  (8) 

j=l  3  3 


where 


Hjfo) 


Pj  *=  probability  of  a  j  cloud  type 
=  arrival  rate  of  cloud  type  j 


F[hj (t) ] 


Figure  11  shows  the  PSD  of  the  combination  of  cloud  sizes  from  200  m  to  5  km 
moving  at  a  constant  velocity  of  10  m/sec  across  a  200  m  pixel.  The  result¬ 
ing  spectra  as  expected  is  intermediate  between  those  shown  in  Figures  5-9 
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4 .  CONCLUSIONS 

The  results  of  these  calculations  show  that  moving  clouds  under 
the  conditions  described  in  this  paper  represent  a  rather  severe  clutter 
problem  for  LWIR  downward  systems.  For  system  footprints  on  the  order 
of  200  m  the  clu-ter  noise  would  vary  over  the  range  4  kw/sr  v*®  to 
40  kw/sr  ym. 

The  work  performed  above  represents  a  first  order  approximation 
to  cloud  clutter  at  LWIR  wavelengths.  For  a  more  refined  model  one  must 
consider  realistic  cloud  situations.  Instead  of  taking  an  average  of 
cloud  sizes  as  above,  one  should  more  properly  expect  the  occurrence  of 
(at  lease  for  cumulus  clouds)  larger  clouds  to  decrease  exponentially.^ 
Still  another  modification  is  needed  to  correct  for  the  spatial  transfer 
function  of  the  sensor.  This  would,  of  course,  have  more  of  an  effect 
on  smaller  sizes. 


^Palnk,  V.G.,  The  Size  Distribution  of  Cumulus  Clouds  in  Representative 
Florida  Populations,  J.  Appl.  Meteorology,  Vol.  8,  1969. 
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FIGURE  1 

Hypothetical  radiance  trace  observed 
by  a  down  looking  sensor  for  a  standard 
atmosphere  with  clouds  at  altitude  of  3  and  6  km. 
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A  typical  sample  of  Poisson 
Arrival. 
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FIGURE  3 

Assumed  radiance  trace  waveforms 
for  the  Paper 
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The  PSD  of  10  laa  clouds  moving  at  30  m/sec  across  a  1  km 
pixel  using  waveform  Fig.  3C. 
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The  PSD  of  combining  cloud  sizes  from  200  m 
co  5  tan  moving  at  10  m/ser  across  a  200  m 
pixel  using  waveform  Fig.  3B. 


AFGL  MAGNETOMETER  NETWORK 


A  Program  Overview 


The  goal  of  this  task  is  to  use  the  AFGL  Magnetometer  Network 
data  base  to  improve  the  Air  Force's  ability  to  specify  and  predict 
geomagnetic  activity.  This  information  is  very  useful  in  the  studying 
of  propagation  of  electro-magnetic  waves  and  messages  in  the  atmosphere. 
Initial  efforts  are  aimed  at  generating  an  easily  accessible  data  base 
for  1-2  years  of  data.  Subsequent  work  will  use  the  data  in  spectral 
and  maximum  entropy  modelling. 

Projected  users  of  this  data  and  analysis  include  Headquarters  of 
Air  Weather  Service,  AFSATCOM,  and  SAC. 


30 


INTRODUCTION 


Space  affects  operational  systems  because  it  is  a  dynamic  environment 
with  daily  and  seasonal  variations  and  with  naturally  occurring  disturbances. 
A  complete  understanding  of  naturally  occurring  disturbances  in  the  magneto¬ 
sphere  and  the  ionosphere  requires  an  understanding  of  the  solar  contribution 
to  such  disturbances.  The  Air  Force  is  concerned  with  the  degrading  effects 
on  its  communications  and  surveillance  systems  which  result  from  interaction 
of  those  systems  with  a  dynamic  space  environment.  However,  the  occurrence 
of  the  magnetic  storms  or  proton  showers  which  cause  these  disturbances 
cannot  be  predicted  without  a  knowledge  of  the  solar  activity  which  causes 
them. 


During  periods  of  high  solar  activity,  a  greatly  enhanced  charged 
particle  population  propagates  from  the  sun  through  interplanetary  medium. 

The  solar  energetic  particles  emitted  from  the  sun  on  interplanetary  magnetic 
field  lines  leading  from  the  sun  to  the  vicinity  of  the  earth  will  travel 
through  the  earth's  magnetosphere  and  impinge  on  the  polar  ionosphere.  Solar 
particle  events  have  a  deleterious  effect  not  only  on  polar  communication 
systems  but  also  on  satellite  sensors  that  are  irradiated  by  solar  particle 
fluxes . 


AEROSPACE  MAGNETIC  MONITORING 


The  ability  to  specify  magnetic  activity  levels  in  real  time  and  to 
predict  such  activity  hours  in  advance  is  urgently  needed  to  support  various 
Air  Force  agencies,  including  AFSATCOM,  SAC,  NORAD  AND  AFGWC.  Other  users 
would  include  SAMSO,  ARPA  and  the  various  groups  engaged  in  fuel  and  mineral 
surveys.  A  program  for  specification  and  prediction  requires  a  system  for 
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the  collection  of  magnetic  activity  data,  its  transmission  to  a  central 
station  and  the  means  for  near  real-time  operation  on  the  data.  This  is 
being  done  by  the  AFGL  magnetic  monitoring  network. 


Geographical  locations  of  the  stations  in  the 
AFGL  Magnetometer  Network. 


Severe  magnetospheric  disturbances,  called  magnetosphere  storms, 
produce  a  variety  of  adverse  effects  on  operational  military  systems.  The 
effects  result  from  severe  departures  from  the  normal  ionospheric  and 
magnetospheric  conditions  upon  which  systems  such  as  military  communications 
depend  for  operation.  The  objectives  of  the  AFGL  Magnetometer  Network 
program  include  the  development  of  methods  to  predict  magnetospheric  storms 
in  order  to  lessen  their  impact  on  military  operational  capability  and  the 
provision  of  real  time  data  to  operational  military  systems  which  require  a 
knowledge  of  magnetospheric  conditions  for  optimum  performance. 

A  magnetospheric  storm  produces  a  number  of  phenomena;  its  manifesta¬ 
tion  as  a  disturbance  of  uhe  geomagnetic  field  is  called  a  magnetic  storm. 
The  AFGL  program  chooses  the  magnetic  field  as  the  parameter  to  be  studied 
because  of  its  direct  dependence  on  the  basic  processes  producing  the 
disturbance  and  the  relative  ease  with  which  it  may  be  measured  by  ground 


magnetometers.  The  practical  approach  chosen  is  to  measure  magnetic- 
field  disturbances  at  a  number  of  selected  locations  in  the  United  States 
with  magnetometers  which  cover  an  appropriate  range  of  amplitude  and 
frequency  and  to  transmit  the  resulting  data  in  real  time  to  a  central 
location.  The  data  obtained  are  measurements  of  the  vector  magnetic  field 
at  intervals  of  one  second  and  measurements  of  the  vector  time  derivative 
of  the  field  for  fluctuations  with  periods  between  one  and  1000  seconds. 

Each  instrumented  data-collection  station  (DCS)  operates  continuous¬ 
ly  and  automatically,  unattended  except  for  routine  maintenance.  Five 
stations  are  located  approximately  along  the  line  of  55°N  corrected  geo¬ 
magnetic  latitude,  and  two  additional  stations  lie  along  the  line  of  40°N 
corrected  geomagnetic  latitude.  Data  from  each  station  are  returned  in 
real  time  on  commercial  voice-grade  communication  circuits  to  a  single 
data-acquisition  station  (DAS)  located  at  the  Air  Force  Geophysics  Labora¬ 
tory  at  Hanscom  Air  Force  Base,  Massachusetts.  The  DAS  performs  real  time 
processing,  reduction,  and  display  of  the  data  and  stores  processed  data 
in  a  permanent  file  for  subsequent  analysis.  All  facilities  involved  are 
dedicated  to  the  program,  so  essentially  uninterrupted  operation  over  an 
extended  time  period  is  possible. 

The  important  features  of  the  magnetometer  network  are  (a)  the 
ordered  array  of  stations,  (b)  the  measurement  of  both  the  vector  field  and 
its  time  derivative  at  a  rapid  sampling  rate  by  identical  instruments  pro¬ 
ducing  directly  comparable  data,  (c)  the  transmission  of  these  data  in 
real  time  to  AFGL,  (d)  the  automatic  real  time  processing. 

All  data  words  are  11  bits  in  length,  with  four  parity  bits  added 
for  error— correction  purposes.  Thus,  a  frame  of  data,  as  prepared  for 
transmission,  consists  of  240  15-bit  words;  these  are  stored  in  a  matrix 
constructed  of  15  240-bit  shift  registers.  Since  one  sampling  interval 


follows  the  preceding  one  without  pause,  while  data  from  each  interval 
need  to  be  stored  until  the  assigned  transmission  time  (between  zero  and 
ten  seconds  later),  two  identical  matrices  are  used  in  alternation,  each 
being  filled  in  one  ten-second  interval  and  storing/transmitting  in  the  next. 
The  output  is  a  serial  stream  of  3600  bits  (240  words  of  15  bits  each) 
clocked  out  at  a  rate  of  4800  bits/sec;  the  duration  of  the  bit  stream  is 
therefore  750  ms. 

In  addition  to  the  redundant  coding  obtained  by  using  four  parity 
bits,  a  technique  of  interleaving  bits  from  all  the  words  is  used  to  reduce 
the  effect  of  a  noise  burst  lasting  longer  than  the  duration  of  one  bit. 

The  first  bit  of  each  word  in  sequence  is  transmitted,  then  the  second  bit 
of  each  word,  and  so  on  until  the  frame  is  complete.  Bursts  of  errors  are 
thus  spread  out  to  cause  no  more  than  one  bit  error  per  word  (which  is  fully 
correctable)  unless  a  burst  persists  for  longer  than  50  ms.  Errors  in  the 
data  from  line  noise  are  entirely  negligible  in  practice. 


THE  COMMUNICATION  LINK 

The  data-acquisition  station  at  AFGL  is  linked  to  the  data-collection 
stations  by  a  voice-grade  communications  network  leased  from  Western  Union. 
The  network  is  full-duplex  (i.e.  affording  simultaneous  two-way  communica¬ 
tion,  an  outlink  and  an  inlink).  It  consists  of  nine  segments  laid  out  to 
form  a  tree,  as  shown  with  AFGL  at  the  base  and  a  DCS  at  the  end  of  each 
branch.  An  outbound  signal,  therefore,  can  be  sent  simultaneously  to  all  of 
the  DCS's,  but  they  must  respond  in  turn.  At  each  hub  the  first  branch  with 
traffic  seizes  the  line  and  holds  it  until  the  transmission  is  completed. 
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DATA  CONTROL 


The  functions  of  the  DAS  may  be  grouped  in  two  areas:  network  control, 
discussed  in  this  section,  and  data  reception  and  processing.  A  schematic 
representation  of  the  DAS  is  shown.  Ten  stations  may  be  accommodated  by 
the  communications  system.  The  inlink  is  used  for  data  return,  time  shared 
by  all  of  the  DCS's,  each  of  which  transmits  a  frame  of  digitized  data  in 
a  programmed  sequence.  Each  DCS  responds  according  to  instructions  contain¬ 
ed  in  the  network-control  signal,  which  synchronizes  the  taking  of  data 
samples  by  the  scientific  instruments  and  the  transmission  of  these  data  at 
the  proper  time.  Signal  propagation  delays  are  precisely  compensated  so 
that  sampling  times  at  all  stations  are  simultaneous  to  within  several 
milliseconds.  The  details  of  sampling,  storing,  and  transmitting  follow  a 
stored  microprocessor  program  of  the  DCS. 

REPRESENTATION  OF  DAS 
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Until  the  past  year,  programs  which  used  the  data  needed  to  incorporate 
a  means  of  checking  and  editing  the  raw  data.  This  proved  to  be  a  very 
difficult  and  time  consuming  activity. 

Bedford  Research  Associates  has  written  a  series  of  computer  routines 
which  successfully  unpack  this  very  complex  data  base.  Now,  for  the  first 
time,  AFGL  scientists  have  access  to  processed  data  with  relative  eash  (with 
contractor  support)  and  meaningful  investigations  are  now  beginning.  The 
project  is  at  a  critical  phase.  Any  interruption  of  support  or  change  in 
personnel  would  cause  substantial  delays  for  AFGL  scientists.  While  the 
system  is  well  documented,  it  is  also  one  that  will  require  many  months  of 
time  and  effort  before  someone  else  could  be  at  a  stage  to  eveh  think  about 
using  the  programs. 

This  processing  capability  will  be  used  to  generate  basic  user  tapes. 
These  tapes,  containing  edited  data,  processed  data,  will  be  arranged  in 
compressed  blocks  covering  one  minute  intervals.  The  basic  user  tapes  will 
be  made  available  to  the  scientific  community. 


REALTIME  DATA 

The  network  minicomputer  has  the  capability  of  doing  some  data 
reduction  in  real  (or  near-real)  time  and  automatically  providing  the  results 
to  users  via  an  appropriate  communication  link.  The  first  such  system  is 
under  construction;  it  will  provide  special  15-rainute  alerts  and  90-minute 
magnetometer  ranges  from  all  stations  to  the  Air  Weather  Service  (AWS)  at 
the  Offutt  AFB  near  Omaha,  NB,  via  an  existing  teletype  network.  In  this 
system,  the  minicomputer  delivers  data  every  15  minutes  to  a  fully  automated 
microprocessor-based  responder,  which  constructs  formatted  messages  and 
transmits  them  to  the  AWS  computer.  It  is  expected  that  a  realtime  detection 
of  sudden  commencements  and  impulses  will  also  be  incorporated  into  this 
system. 
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DATA  ACCESSING  ROUTINE 


Since  a  large  data  base  is  being  created,  a  standard  routine  was 
envisaged  to  access  any  point  of  the  data.  This  would  eliminate  any 
redundant  efforts  needed  to  access  separate  parts  of  the  data  in  separate 
analyses.  However,  detailed  consideration  was  necessary  to  provide  a  sub¬ 
routine  which  would  be  easily  usable  in  applications. 

The  data  are  read  from  file  MSDATA  by  functions  MGREAD  (MSINN,  ID, 
IDATA,  NDATA,  NSTA) .  MSINN  is  the  earliest  acceptable  time  in  milliseconds 
since  1970.  If  positive,  the  next  satisfactory  data  set  is  provided.  If 
MSINN  is  -1,  the  last  data  set  which  was  extracted  is  reused.  If  a  milli¬ 
second  count  is  set  negative,  then  only  that  number  of  milliseconds  or  the 
closest  time  after  that  specified  is  acceptable.  If  only  a  particular 
station  is  wanted,  the  negative  of  the  station  number  is  entered  in  NSTA: 
positive  values  are  not  checked. 

The  data  set  is  specified  by  MSINN  and  possibly  NSTA.  ID  indicates 
which  data  are  to  be  extracted  from  the  set. 

Of  the  two  decimal  digits  of  ID,  the  first  specifies  the  type  of  data, 
the  second  subdivides  the  type.  The  values  34  and  14  give  the  vector  values 
of  the  magnetic  field  each  second  and  the  change  in  the  magnetic  field  each 
.2  seconds. 

The  general  requirements  defining  the  main  program  were  flexibility  in 
use,  the  possibility  of  expansion  and  modification  and  the  capability  to  be 
used  in  either  batch  mode  of  interactively. 

Flexibilities  can  be  obtained  through  use  of  the  modular  structure; 
each  input  card  specifies  a  modeule  to  be  executed.  This  allows  modules  to 
be  added  or  modified  individually  and  puts  the  ordering  of  modular  operations 
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under  control  of  the  input  deck.  For  ease  in  interactive  use  and  in  batch, 
the  input  is  free  form  with  a  keyword  being  a  mnemonic  in  the  operation  to 
be  done.  Interactive  use  restricts  the  size  of  the  allowed  program  requiring 
an  overlay  structure.  The  main  overlay  was  defined  to  be  small  as  possible 
with  a  higher  overlay  to  interpret  the  keyword  cards  and  control  the  main 
overlay.  This  allows  the  other  overlays  to  be  treated  as  modules.  The 
separate  overlays  act  as  individual  programs  except  that  (1)  certain  data 
can  be  passed  between  them  or  saved  and  reused  and  (2)  they  can  link  each 
other  without  requiring  another  keyword  card. 

The  data  is  packed  into  12  bit  quantities,  which  will  be  called  slots 
to  avoid  conceptual  problems  with  16  and/or  12  bit  words.  The  slots  are 
grouped  into  sets;  the  sets  in  turn  are  grouped  into  records. 

The  data  stored  are  signed  integers  in  2's  compliment  form.  A  slot 
may  contain  a  single  11-bit  integer,  sign  and  a  10-bit  magnitude,  plus  a 
flag  in  the  low  order  bit,  which  is  set  if  the  data  is  suspect.  The  initial 
slots  in  a  set,  which  contain  information  pertinent  to  accessing  the  data, 
are  12-bit  signed  integers;  i.e.,  they  do  not  have  a  flag  as  the  low  order 
bit.  However,  a  string  of  data  may  be  compressed  into  an  initial  value,  a 
string  of  deltas  (the  increment  from  the  present  value  to  the  next),  and  the 
final  value.  The  compression  results  from  the  restriction  of  the  deltas  to 
6, A, 3, 2,  or  1/bit  quantities  which  are  packed  into  the  slots  giving  packing 
factors  of  2, 3, A, 6  and  12.  The  shortened  deltas  are  packed  from  low  order 
to  high  which  is  opposite  many  conventions.  The  first  delta  is  inaccurate. 

To  reconstruct  the  data  the  deltas  must  be  used  in  reverse  order  to  produce 
the  data  from  the  next  to  the  last  value  to  the  second.  In  the  case  of  a 
value  remaining  constant,  no  deltas  are  stored  and  only  one  value  is  stored. 

Some  slots  may  contain  strings  of  one-bit  flags.  In  this  case  the 
packing  is  from  high  order  to  low. 
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The  first  three  slots  of  a  data  set  gives  the  type  of  data  set,  the 
subtype  of  data  set,  and  the  total  number  of  slots  comprising  the  data  set. 

A  data  set  has  two  different  formats  depending  on  whether  packing  was 
used.  The  fourth  word  is  set  to  slots  containing  the  information  needed  to 
access  the  data  in  the  two  cases.  Brackets,  [],  indicate  data  which  might 
be  packed.  Vector  data  is  given  in  three  arrays:  all  x-data  first,  then  y, 
then  z.  X  is  north,  y  is  east,  and  z  is  down. 

Since  the  records  have  a  maximum  length  of  2560  slots  and  the  sets  of 
slots  have  a  variable  length,  a  data  type  of  -1  terminates  the  records. 

There  may  be  slots  filled  with  garbage  after  the  -1  data  type  but  no  data 
should  follow  until  the  next  record. 

The  tape  ends  with  a  double  end  of  file,  a  single  end  of  file,  or  the 
physical  end  of  the  tape.  Since  some  tapes  have  been  restarted,  an  end  of 
file  may  be  followed  by  subsequent  data. 

In  addition  to  the  transmission  errors  which  are  detected  and  either 
corrected  or  flagged,  errors  of  various  types  evidence  themselves  in  the  time 
parameter;  among  these  are:  erroneous  time,  data  apparently  out  of  chrono¬ 
logical  sequence,  and  subsequent  data  overwriting  the  start  of  a  data  tape. 


PROGRAM  STRUCTURE 

The  program  was  designed  for  an  overlay  structure  to  save  space.  The 
main  overlay  was  designed  as  simply  as  possible;  its  sole  function  is  to  call 
higher  overlays  upon  command  and  pass  data.  Significant  data  is  passed 
via  commons  as  are  the  commands  controlling  the  calling  of  overlays. 
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A  few  generally  useful  subroutines  are  included  which  are  used  by  several 
overlays.  One  results  of  this  organization  is  that  general  data  read  to 
control  one  keyword  may  be  available  for  others  even  when  different  overlays 
are  used. 

The  program  was  constructed  and  checked  under  a  system  which  kept  a 
record  of  all  modifications. 

Overlay  (1,0) - Control 

Overlay  (1,0)  controls  the  program.  Keyword  cards  are  read  in  and  act¬ 
ed  upon.  For  ease  in  use  these  are  free  field  cards  modelled  after  the  CDC 
job  control  cards,  except  that  numeric  fields  may  be  either  integer  or  real. 
Keywords  which  change  parameters  or  do  simple  tasks  are  acted  upon  in  this 
overlay.  Thos  needing  magnetogram  data  request  the  main  overlay  to  call  in 
the  appropriate  overlays  before  returning  control  to  this  overlay.  In  this 
way,  the  modular  structure  of  the  entire  program  is  maintained;  control  of 
the  sequence  of  operations  is  given  to  the  person  running  the  program  through 
the  keyword  cards. 

Overlay  (2,0) - Data  Base  Access 

In  many  ways  this  overlay  does  the  main  task  of  the  program:  a  data 
tape  is  read  and  a  compact  data  file  is  created.  The  data  tape,  MGDATA,  is 
read  by  MGREAD  as  described  in  section  2.1.1.  The  data  is  checked  for  bad 
spots  while  being  unpacked  and  bad  data  is  replaced  by  99990.  or  99999. 

(This  could  be  changed  through  common  /MGREADC/) .  The  data  times  returned 
are  checked  if  the  data  is  to  be  used.  Times  showing  jumps  greater  than  15 
seconds  are  stored  until  either  the  time  discrepancy  resolves  itself  or  the 
error  storage  area  is  filled.  In  the  former  case,  the  times  of  the  saved 
data  is  corrected  before  the  data  is  used;  in  the  latter,  the  new  time  is 
considered  correct.  Time  reversals  due  to  long  sections  of  bad  data  impede 
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the  use  >  1  owing  data.  For  this  reason  MGDATA  is  not  normally  rewound 

before  use  Si.<  that  the  tape  may  be  positioned  externally  from  the  reading 
routine . 

Once  the  data  has  been  unpacked  and  checked  for  problems,  especially 
time  discrepancies,  and  it  is  from  one  of  the  stations  and  time  interval 
desired,  the  data  is  stored  by  routine  STODTA  in  the  buffer  1BUF.  When 
this  buffer  is  full,  the  routine  AVGOUT  is  called  to  pack  the  data  into  the 
buffer  I0B.  The  data  is  packed  by  averaging  controlled  by  two  main  para¬ 
meters:  NMSAVG  and  MSDAVG.  NMSAVG  is  the  number  of  milliseconds  to  be 
averaged  over.  MSDAVG  is  the  displacement  from  one  average  to  the  next. 

NTAVG  AND  NTDAVG  are  the  number  of  time  points  in  the  average  and  in  the  ' 
displacement  assuming  MSPT  milliseconds  per  time  point.  Depending  on  the 
relative  values  of  these  parameters,  the  averages  may  overlap  producing 
smoothing  or  they  may  be  widely  separated  for  sampling.  Averaging  over 
one  sample  results  in  direct  use  of  the  data.  Data  which  was  flagged  by 
MGREAD  is  not  used  and  averages  with  no  data  are  flagged  by  AVGFLG,  set  to 
65536.  Times  having  no  data  are  flagged  by  FFLAG,  also  set  to  65536.  Data 
not  needed  for  further  averaging  is  flushed  from  the  buffer  and  the  remaining 
data  repositioned. 

When  the  output  buffer  is  filled,  subroutine  I0BUT  is  called  to  output 
the  data  as  a  record  of  file  AVGEN  and  to  adjust  the  buffer  parameters  to 
allow  more  data  to  be  placed  in  it. 


Overlay  (3,0) - Plotting 

This  overlay  extracts  data  from  AVGEN  and  organizes  it  for  plotting. 

The  plotting  request  can  be  for  any  length  of  time,  and  for  any  number  of 
stations.  To  allow  for  this  flexibility,  the  routine  may  break  a  requested 
time  period  into  several  plots.  Each  plot  is  limited  to  MNP  points,  current¬ 
ly  equal  to  720.  This  may  result  in  sampling  the  points  rather  than  use  all 
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of  them;  in  this  routine,  sampling  rather  than  averaging  is  done. 

This  overlay  and  the  next  can  be  used  interactively  to  examine 
sections  of  the  data  in  as  much  detail  as  desired. 


Overlay  (4,0) - Specialized  Storage 

This  overlay  transfers  data  between  the  standard  file  and  a 
random  access  main  storage  file.  When  data  is  transferred  from  AVGEN, 
the  two  files  do  not  have  to  agree  in  their  parameters  controlling  the 
spacing  of  data  points.  Data  from  AVGEN,  falling  with  the  time  interval 
and  from  the  desired  stations  will  be  transferred  to  the  position  in 
MSFILE  whose  time  is  the  closest.  When  a  new  AVGEN  file  is  created,  its 
time  parameters  are  derived  from  thos  of  MSFILE.  To  make  the  file  AVGEN 
more  compact,  only  data  from  the  specified  stations  are  used  and  only 
times  having  at  least  one  valid  piece  of  data  are  used. 
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MAGNETOGRAM  PROCESSING 


0.  Introduction 

The  magnetogram  data  were  collected  at  seven  stations,  transmitted 
to  the  Air  Force  Geophysical  Laboratory,  and  stored  on  magnetic  tape.  The 
data  from  one  station  consists  of  readings  of  the  magnetic  field  of  the 
earth,  the  rate  of  change  of  the  field,  and  some  engineering  data.  This 
data  is  collected,  packed,  and  stored  for  10  second  intervals  at  each  site. 
The  packed  data  is  transmitted  in  ten  second  blocks  from  the  station  to 
AFGL.  A  dedicated  mini-computer  reorganizes  the  data  and  stores  it  on 
magnetic  tapes.  These  tapes  form  the  data  base. 

The  stations  are  listed  in  Table  1.  More  detailed  information  on  the 
magnetometer  data  network  may  be  found  in  D.J.  Knecht,  R.O.  Hutchinson,  and 
C.W.  Tsacoyeanes,  AN  INTRODUCTION  TO  THE  AFGL  MAGNETOMETER  NETWORK,  Air 
Force  Geophysics  Laboratory,  April,  1979.  (Unpublished). 

The  main  ideas  governing  the  construction  of  a  routine  to  read  this 
data  base  and  of  a  main  program  to  analyze  and  display  the  data  are  given 
in  section  1.  Section  2  covers  aspects  of  the  analysis  program  as  well  as 
detailed  descriptions  of  the  data  bases  used  or  created;  in  particular  the 
structure  of  the  main  data  base  is  given  in  section  2.1.1.  Particulars  of 
the  use  of  the  program  are  in  section  3. 

1.  Basic  Concepts 

1.1  Data  Accessing  Routine 

Since  a  large  data  base  is  being  created,  a  standard  routine  was 
envisaged  to  access  any  point  of  the  data.  This  would  eliminate  any 
redundant  efforts  needed  to  access  separate  parts  of  the  data  in  separate 
analyses.  However,  detailed  consideration  was  necessary  to  provide  a  sub¬ 
routine  which  would  be  easily  usable  in  applications. 
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The  data  are  read  from  file  MSDATA  by  function  MGREAD,  (MSINN,  ID, 
IDATA,  NDATA,  NSTA) .  MSINN  is  the  earliest  acceptable  time  in  milliseconds 
since  1970,  If  positive,  the  next  satisfactory  data  set  is  provided.  If 
MSINN  is  -1,  the  last  data  set  which  was  extracted  is  reused.  If  a  milli¬ 
second  count  is  set  negative,  then  only  that  number  of  milliseconds  or  the 
closest  time  after  that  specified  is  acceptable.  If  only  a  particular 
station  is  wanted,  the  negative  of  the  station  number  is  entered  in  NSTA; 
positive  values  are  not  checked. 

The  data  set  is  specified  by  MSINN  and  possibly  NSTA.  ID  indicates 
which  data  are  to  be  extracted  from  the  set. 

Of  the  two  decimal  digits  of  ID,  the  first  specifies  the  type  of 
data,  the  second  subdivides  the  type.  The  possibilities  are  listed  in 
table  1.2.  The  values  34  and  14  give  the  vector  values  of  the  magnetic 
field  each  second  and  the  change  in  the  magnetic  field  each  .2  seconds. 

Additional  controls  are  provided  through  commons.  In  particular, 
/MGREDC/  and  /STANUM/  provide  switches  for  optional  handling  of  bad  data; 
/MGDIMF/  indicates  the  present  location  on  the  tape  and  control  restarting; 
/MGDATC/  control  the  unit  read. 

1.2  Program  Requirements 

The  general  requirements  defining  the  main  program  were  flexibility 
in  use,  the  possibility  of  expansion  and  modification  and  the  capability 
to  be  used  in  either  batch  mode  or  interactively. 

Flexibility  can  be  obtained  through  a  modular  structure  where  each 
input  card  specifies  a  module  to  be  executed.  This  allows  modules  to  be 
added  or  modified  individually  and  puts  the  ordering  of  modular  operations 
under  control  of  the  input  deck.  For  ease  in  interactive  use  and  in  batch, 
the  input  is  free  form  with  a  keyword  being  a  mnemonic  for  the  operation  to 
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be  done.  The  interactive  use  restricts  the  size  of  the  allowed  program 
requiring  an  overlay  structure.  The  main  overlay  was  to  be  as  small  as 
possible  with  a  higher  overlay  to  interpret  the  keyword  cards  and  to 
control  the  main  overlay.  This  allows  the  other  overlays  to  be  treated 
as  modules.  The  separate  overlays  act  as  individual  programs  except  that 
certain  data  can  be  passed  between  them  or  saved  and  reused  and  that  they 
can  link  to  each  other  without  requiring  another  keyword  card. 


45 


a 

o 

'H 

W 

CO 


«  X 
ctf  C  <V 
Vi  <0  *0 

(0  u  c 

OHM 


m  n  »  co  N 


<D 

"O  -X) 
*-»  P 
co  m 


ID 
r— 4 

tfl 

•H 

Q 


U0  NO 


<  O  M  M 
3  co  26  X 


3  * 


u 


CO 

2 

O 

M 

H 

< 

H 

co 

Pi 

W 

H 

S 


cfl 


< 

Q 

M 

H 

<C 

.J 

< 

o 

H 

3 

CO 

3 

fc 

* 

u* 

O 

H 

•H 

W 

•U 

*» 

•* 

•» 

*» 

•» 

«K 

* 

Z 

t-4 

3 

o 

CO 

-J 

PQ 

c 

o 

O 

a 

w 

a 

a 

u 

Pl, 

o 

< 

V-4 

2 

u 

s: 

CO 

H 

►J 

£ 

oc 

c 

vD 

co 

cm 

oo 

ON 

ON 

sO 

o 

• 

• 

• 

• 

• 

• 

• 

a\ 

a* 

<3 

i — 4 

vy 

o 

On 

t— 4 

ro 

<r 

o 

W 

cm 

ro 

ro 

m 

ro 

ro 

u 

•H 

4-> 

0)  4-J 

P  <0 

CM 

H 

ro 

00 

00 

a* 

CM 

6C 

• 

• 

• 

• 

• 

• 

• 

cfl 

in 

<r 

NO 

m 

m 

O 

O 

S  z 

m 

in 

m 

m 

m 

-3“ 

00 


c 

* — 4 

H 

ro 

ON 

ro 

m 

sO 

r— 4 

o 

• 

• 

• 

• 

• 

• 

• 

03 

►J 

A* 

ro 

O 

eg 

r — 4 

CM 

o 

O 

H 

O 

Os 

00 

a» 

00 

CM 

v-4 

3 

f— 4 

* — ) 

i-H 

c 

x: 

O 

D. 

-H 

0J 

Vi 

W 

Vi 

00 

03 

ro 

CM 

o 

NO 

CM 

00 

A 

U 

o 

►J 

• 

• 

• 

• 

• 

• 

• 

o 

03 

co 

<r 

<3- 

CM 

CM 

A. 

kJ 

o 

z 

<r 

-<r 

<1- 

<3“ 

CO 

46 


ID  Data  Item  Total  Number 

Value  Every  sec.  of  Values 


rH 

rH 

rH 

H 

X 

X 

X 

X 

cO 

*H 

*H 

•H 

•H 

X 

X 

X 

X 

3 

3 

3 

3 

co 

44 

o 

O 

O 

O 

X 

X 

X 

X 

3 

3 

3 

3 

rH 

rH 

rH 

rH 

4J 

CO 

u 

U 

O 

O 

3 

3 

3 

3 

rH 

H 

rH 

rH 

Eh 

Eh 

Eh 

Eh 

CO 

P 

rH 

*H 

H 

rH 

u. 

Eh 

Eh 

P 

42 

42 

42 

42 

Eh 

Eh 

f* 

Eh 

0) 

a> 

a 

V 

60 

U 

U 

O 

O 

rH 

rH 

rH 

»H 

CO 

CO 

CO 

CO 

60 

c 

U 

M 

M 

0) 

0) 

0) 

0) 

CO 

cO 

(0 

tO 

P 

P 

p 

p 

c 

•H 

<0 

co 

CO 

CO 

c 

C 

3 

C 

44 

44 

44 

4J 

cO 

CO 

C0 

CO 

«H 

u 

<u 

QJ 

<U 

0) 

■H 

•H 

*H 

■H 

o 

O 

o 

O 

o 

o 

o 

O 

c0 

p 

Q> 

co 

CO 

co 

CO 

Eh 

Eh 

Eh 

Eh 

H 

H 

H 

H 

V 

o 

o 

o 

44 

0) 

<U 

co 

0) 

c 

P 

d 

•H 

4J 

•M 

44 

M 

44 

•u 

44 

u 

44 

44 

44 

P 

44 

44 

44 

p 

*H 

60 

c 

c 

c 

o 

c 

C 

c 

o 

c 

c 

c 

O 

c 

(2 

c 

o 

rH 

(40 

d 

Q) 

0) 

0) 

44 

<u 

<u 

0) 

4-1 

0) 

0) 

<u 

44 

0) 

cu 

a> 

44 

CO 

c 

E*2 

c 

c 

c 

u 

c 

c 

C 

o 

c 

c 

3 

o 

c 

e 

(2 

o 

C 

w 

o 

o 

o 

a) 

o 

o 

O 

0) 

o 

o 

0 

a) 

o 

o 

O 

(1) 

o 

rH 

CL 

CL 

Cl 

> 

9* 

CL 

3.  > 

cl 

CL 

CL  > 

CL 

a 

CL  > 

»H 

60 

CO 

B 

e 

a 

e 

a 

a 

e 

6 

e 

6 

e 

E 

44 

O 

44 

O 

o 

o 

rH 

o 

o 

o 

iH 

o 

o 

o 

iH 

o 

o 

O 

rH 

*H 

rH 

*H 

a 

o 

u 

rH 

o 

o 

O 

rH 

o 

u 

u 

rH 

u 

u 

<4 

rH 

TJ 

«0 

60 

3 

3 

3 

3 

T3 

d 

*H 

X 

N 

Pm 

X 

>s 

N 

Eh 

X 

Jo 

N 

Eh 

X 

Jo 

N 

Eh 

< 

-< 

P 

rH 

CM 

<o 

rH 

CM 

CO 

<r 

iH 

CM 

CO 

<r 

rH 

CM 

CO 

■sf 

rH 

CM 

H 

«H 

rH 

CM 

CM 

CM 

CM 

co 

CO 

CO 

CO 

<* 

m 

m 

47 


i 

A 


DATA  TO  BE  EXTRACTED 


2. 


Struc  tures 


The  main  structures  are  given  in  this  section 

2.1  Data  Structures 

2.1.1  Magnetogram  Data  Tapes 

The  standard  file  name  for  these  is  MGDATA.  The  tapes  are  7  track, 

800  b.p.i.  A  short  description  of  the  format  is  given  belovr;  for  a  complete 
description,  see  D.J.  Knecht,  Archive  Data  Tapes  of  the  AFGL  Magnetometer 
Network,  Air  Force  Geophysical  Laboratory,  to  be  printed  (Unpublished). 

The  data  are  packed  into  12  bit  quantities,  which  will  be  called  slots 
to  avoid  conceptual  problems  with  16  and/or  12  bit  words.  The  slots  are 
grouped  into  sets;  the  sets  in  turn  are  grouped  into  records. 

The  data  stored  are  signed  integers  in  2's  compliment  form.  A  slot 
may  contain  a  single  11-bit  integer,  sign  and  10-bit  magnitude,  plus  a  flag 
in  the  low  order  bit,  which  is  set  if  the  data  is  suspect.  The  initial 
slots  in  a  set  which  contain  information  pertinent  to  accessing  the  data, 
are  12-bit  signed  integers;  i.e.  they  do  not  have  a  flag  as  the  low  order  bit. 
However,  a  string  of  data  may  be  compressed  into  an  initial  value,  a  string 
of  deltas  (change  from  the  present  value  to  the  next),  and  the  final  value. 

The  compression  results  from  the  restriction  of  the  deltas  to  6, 4, 3, 2,  or  1 
bit  quantities  which  are  packed,  the  slots  giving  packing  factors  of  2, 3, 4, 6, 
and  12.  The  shortened  deltas  are  packed  from  low  order  to  high,  which  is 
opposite  many  conventions.  The  first  delta  is  inaccurate.  To  reconstruct  the 
data  the  deltas  must  be  used  in  reverse  order  to  produce  the  data  from  the 
next  to  the  last  value  to  the  second.  In  the  case  of  a  value  remaining 
constant,  no  deltas  are  stored  and  only  one  value  is  stored. 
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Some  slots  may  contain  strings  of  one-bit  flags.  In  this  case  the 
packing  is  from  high  order  to  low. 

The  first  three  slots  of  a  data  set  define  the  type  of  data  set,  the 
subtype  of  data  set,  and  the  total  number  of  slots  comprising  the  data  set. 

A  data  set  has  two  different  formats  depending  on  whether  packing  was 
used.  The  ordering  of  the  data  is  described  in  Table  2.1.  The  fourth  word 
is  set  to  1  if  packing  was  used  and  to  -1  if  it  was  not.  The  two  columns 
of  slot  numbers  in  Table  2.1  indicate  which  slots  contain  the  information 
needed  to  access  the  data  in  these  two  cases.  Brackets  [],  indicate  data 
which  might  be  packed.  Vector  data  is  given  in  three  arrays:  all  x  data 
first,  then  y,  then  z.  X  is  north,  y  is  east,  and  z  is  down. 

Since  the  records  have  a  maximum  length  of  2560  slots  and  the  sets  of 
slots  have  a  variable  length,  a  daya  type  of  -1  terminates  the  record.  There 
may  be  slots  filled  with  garbage  after  the  -1  data  type  but  no  data  should 
follow  until  the  next  record. 

The  tape  ends  with  a  double  end  of  file,  a  single  end  of  file,  or  the 
physical  end  of  the  tape.  Since  some  tapes  have  been  restarted,  an  end  of 
file  may  be  followed  by  subsequent  data. 

In  addition  to  the  transmission  errors  which  are  detected  and  either 
corrected  or  flagged,  errors  of  various  types  evidence  themselves  in  the 
time  parameter;  among  these  are:  erroneous  time,  data  apparently  out  of 
chronological  sequence,  and  subsequent  data  overwriting  the  start  of  the 
data  tape. 
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TABLE  2 . 1  DATA  FORMATS 


Packing  Used 
Slot 


1 

2 

3 

4 

5,6,7 

8,9.10 

11,12 

13,14 


No  Packing 
Slot 


1 

2 

3 

4 

/ 

/ 

/ 

5,6 


15  7 

16,  on  8,  on 


Quantity 

1  (type) 

0  (subtype) 

Number  of  slots 
Packing  flag:  1  “  packing 
-1  =  no  packing 
No.  of  bits/value  for 
Search  Coil  X,Y,Z 
No.  of  bits/value  for 
Fine  Flux  X.Y.Z 
No.  of  bits/value  for 
Analog  and  Digital  data 
Time  of  receipt  of  data 
in  seconds  past  10  second 
mark  and  milliseconds 
Diagnostic  check  for  mini¬ 
computer 
Actual  data: 

[Search  Coil  X] 

[Search  Coil  Y] 

[Search  Coil  Z] 

[Fine  Flux  X] 

[Fine  Flux  Y] 

[Fine  Flux  Z] 

[Coarse  Flux  X] 

[Coarse  Flux  Y] 

[Coarse  Flux  Z] 


[Analog  data] 
[Digital  data] 

Additional  data 
1 
2 

3 

4 

5 

6,7 

8,9,10 

11,12,13 


50  values  of  each 
Unit=gamma/sec 

10  values  of  each 
Unit=gamma/sec 

3  slots  for  each 
2  values,  unit  =  64 
Flag  string,  align 
with  fine  flux  value 
indicating  location 
of  change 

23  values 
15  values 

13  slots 
5  volt  values 

Temperature  of  electronics 

Temp,  of  sensor 

Station  ID  number 

Serial  Number 

Status  words 

Synch  words 

Time  of  data 
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2.1.2  Averaged  General  Data  File 

The  file  AVGEN  is  the  main  file  for  data  use.  Overlay  (2,0)  creates 
this  file  from  the  data  tapes.  (4,0)  transfers  data  between  this  type  of 
file  and  mass  storage.  (3,0)  plots  data  stored  on  AVGEN.  The  file  OLDAVG 
has  the  same  structure. 

The  first  logical  record  is  a  character  string  labeling  the  file. 

The  first  100  characters  contain  the  date  and  time  which  the  file  was 
created  as  a  standard  file.  Subsequent  sets  of  10  characters  contain  the 
date  and  title  of  each  run  which  modifies  the  file.  A  list  of  abbreviations 
of  the  stations  used  end  the  record. 

The  second  logical  record  consists  of  data  specifying  the  averaging 
and  packing  of  the  magnetic  data  held  on  the  file.  This  consists  of  copies 
of  the  two  commons  /AVGC/  and  /STATC/ .  Each  copy  is  preceeded  by  one  word 
giving  the  length  of  the  common;  this  allows  for  future  changes  without 
invalidating  any  earlier  files.  The  data  in  the  commons  are  described 
individually  in  the  description  of  overlay  (2,0). 

The  following  logical  records  consist  of  as  many  sets  as  can  fit  into 
L10B  words;  LI0B  being  set  to  512.  A  set  consists  of  a  time  word  and  three 
average  data  words  from  each  station  which  is  to  be  included.  Thus,  if  NSTA 
stations  are  used,  a  set  has  3*  NSTA  +  1  words.  The  time  word  is  an  integer 
giving  the  average  time  in  milliseconds  since  1970.  The  data  words  are  the 
average  flux  with  the  value  of  AVGFLG,  set  to  65536.,  indicating  the  lack  of 
data. 


The  data  are  in  chronological  order.  The  number  of  points  per  average 
is  NTAVG  and  the  spacing  of  averages  is  MSDAVG;  thus  the  average  could 
contain  one  point  for  non-averaging,  could  overlap  for  smoothing,  or  could 
be  widely  separated  for  sampling.  The  parameters  may  be  changed  by  input 
cards. 
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2.1.3  Mass  Storage  Data  File 

The  file  MSFILE  is  a  random  access  storage  file  which  allows  large 
amounts  of  data  to  be  stored  and  accessed  easily.  The  data  must  be  put 
into  AVGEN  form  before  using  but  this  can  be  done  quickly  and  without 
imprinting  a  tape. 

The  file  is  a  CDC  random  access  file  with  a  numerical  index.  The 
first  record  is  a  character  string  containing  the  starting  and  ending  date 
and  time,  a  general  label,  a  list  of  station  initials,  and  the  date,  time 
and  title  of  each  run  contributing  to  the  file. 

The  second  record  contains  a  copy  of  common  /RAMXS/ .  The  contents 
are  described  in  detail  in  the  section  on  overlay  (4,0). 

The  other  records  contain  data.  The  first  word  is  the  time  of  the 
first  data  set.  Each  record  contains  data  sets  each  of  which  has  NQSETM 
items:  NDSM  items  from  each  of  NSTAM  stations.  The  time  interval  between 
sets  is  MSDMS  milliseconds.  The  time  of  the  first  data  set  is  given  in 
milliseconds  since  1970.  If  a  datum  is  absent  the  value  FILLMS ,  set  to 
65536.,  is  used.  Each  record  contains  an  interval  of  MSREC  milliseconds 
and  the  initial  record  starts  at  INTMSM  milliseconds  since  1970. 


2.2  Program  Structure 


The  program  was  designed  in  an  overlay  structure  to  save  space. 

The  main  overlay  was  designed  as  simply  as  possible;  its  sole  function 
is  to  call  higher  overlays  upon  command  and  pass  data.  Significant  data 
are  passed  via  commons  as  are  the  commands  controlling  the  calling  of 
overlays.  A  few  generally  useful  subroutes  are  included  which  are  used 
by  several  overlays.  One  result  of  this  organization  is  that  general  data 
read  to  control  one  keyword  may  be  available  for  others  even  when  different 
overlays  are  used. 

The  program  was  constructed  and  checked  under  a  system  which  kept  a 
record  of  all  modifications.  Since  this  system  is  of  more  general  use  it 
is  described  in  appendix  A.l. 

The  overlays  are  described  in  the  following  sections.  The  routines 
are  listed  in  tables  2.2.nR  and  the  commons  in  table  2.2.nC  where  n  is  the 
number  of  the  overlay.  Similarly,  flow  charts  of  the  main  routines  are  in 
figures  2.n. 

2.2.1  Overlay  (1,0)  -  Control 

Overlay  (1,0)  controls  the  program.  Keyword  cards  are  read  in  and 
acted  upon.  For  ease  in  use  these  are  free  field  cards  modeled  after  CDC 
job  control  cards,  except  that  numeric  fields  may  be  either  integer  or  real. 
Keywords  which  change  parameters  or  do  simple  tasks  are  acted  upon  in  this 
overlay.  Those  needing  magnetogram  data  request  the  main  overlay  to  call  in 
the  appropriate  overlays  before  returning  control  to  this  overlay.  In  this 
way,  the  modular  structure  of  the  entire  program  is  maintained;  control  of 
the  sequence  of  operations  is  given  to  the  person  running  the  program  through 
the  keyword  cards. 
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2.2.2  Overlay  (2,0)  Data  Base  Access 


In  many  ways  this  overlay  does  the  main  task  of  the  program:  a  data 
tape  is  read  and  a  compact  data  file  is  created.  The  data  tape,  MGDATA, 
is  read  by  MGREAD  as  described  in  section  2.1.1.  The  data  is  checked  for 
bad  spots  while  being  unpacked  and  bad  data  is  replaced  by  99990.  or  99999. 
(This  could  be  changed  through  common  /MGREDC/) .  The  data  times  returned 
are  checked  if  the  data  are  to  be  used.  Times  showing  jumps  greater  than 
15  seconds  are  stored  until  either  the  time  discrepancy  resolves  itself  or 
the  error  storage  area  is  filled.  In  the  former  case,  the  times  of  the 
saved  data  is  corrected  before  the  data  is  used;  in  the  latter,  the  new  time 
is  considered  correct.  Time  reversals  due  to  long  sections  of  bad  data 
impede  the  use  of  the  following  data;  for  this  reason  MGDATA  is  nut  normally 
rewound  before  use  so  that  the  tape  may  be  positioned  externally  from  the 
reading  routine. 

Once  a  data  set  is  unpacked  and  checked  for  problems,  especially  time 
discrepancies,  and  it  is  from  one  of  the  stations  and  time  interval  desired, 
the  data  is  stored  by  routine  STODTA  in  the  buffer  IBUF.  When  this  buffer 
is  full,  the  routine  AVGOUT  is  called  to  pack  the  data  into  the  buffer  I0B. 
The  data  is  packed  by  averaging  controlled  by  two  main  parameters:  NMSAVG 
and  MSDAVG.  NMSAVG  is  the  number  of  milliseconds  to  be  averaged  over. 

MSDAVG  is  the  displacement  from  one  average  to  the  next.  NTAVG  and  NTDAVG 
are  the  number  of  time  points  in  the  average  and  in  the  displacement 
assuming  MSPT  milliseconds  per  time  point.  Depending  on  the  relative  values 
of  these  parameters,  the  averages  may  overlap  producing  smoothing  or  they 
may  be  widely  separated  for  sampling.  Averaging  over  one  sample  results  in 
direct  use  of  the  data.  Data  which  was  flagged  by  MGREAD  is  not  used  and 
averages  with  no  data  are  flagged  by  AVGFLG,  set  to  65536.  Times  having  no 
data  are  flagged  by  FFLAG,  also  set  to  65536.  Data  not  needed  for  further 
averaging  is  flushed  from  the  buffer  and  the  remaining  data  repositioned. 
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2.2.3  Overlay  (3,0)  Plotting 

This  overlay  extracts  data  from  AVGEN  and  organizes  it  for  plotting. 

The  plotting  request  can  be  for  any  length  of  time,  and  for  any  number  of 
stations.  To  allow  for  this  flexibility,  the  routine  may  break  a  requested 
time  period  into  several  plots.  Each  plot  is  limited  to  MNP  points,  current¬ 
ly  equal  to  720.  This  may  result  in  sampling  the  points  rather  than  use  all 
of  them;  in  this  routine,  sampling  rather  than  averaging  is  done. 

This  overlay  and  the  next  can  be  used  interactively  to  examine 
sections  of  the  data  in  as  much  detail  as  desired. 

2.2.4  Overlay  (4,0)  Specialized  Storage 

This  overlay  transfers  data  between  the  standard  file  and  a  random 
access  main  storage  file.  When  data  is  transferred  from  AVGEN,  the  two 
files  do  not  have  to  agree  in  their  parameters  controlling  the  spacing  of 
data  points.  Data  from  AVGEN,  falling  within  the  time  interval  and  from 
the  desired  stations  will  be  transferred  to  the  position  in  MSFILE  whose 
time  is  the  closest.  When  a  new  AVGEN  file  is  created,  its  time  parameters 
are  derived  from  those  of  MSFILE.  To  make  the  file  AVGEN  more  compact,  only 
data  from  the  specified  stations  are  used  and  only  times  having  at  least 
one  valid  piece  of  data  are  used. 
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START 


FIGURE  2.0 


Overlay  (0,0) 


Base  Overlay 


Overlay  (1,0) 


FICURE  2.1 

Standard  Modular 


Driver/ controller 


Each  overlay  is  one 
or  more  modules. 


FIGURE  2.1.1 


Overlay  (3.0) 


Plotting  Overlay 


FIGURE  2.3 


Overlay  (4,0) 


Conversion  of  File  Types 


(  From  (0,0)  ) 

_ - r~  ^ 

Store  input  quantities 
Initialize  partially 


Complete  initialization 
including  checking  AVGEN 
and  opening  and  checking 
MSFILE 


Set  JI0B  to  point  to 
next  data  set 


(OB:  LIOB  *°  AVGEN ^record^eng^tT 

=rda!Ta-_timel 


READ  next  AVGEN 
record,  adjust 
pointer  JIOB 


Tata  time:  initial  tjjaeZl^ 

Ja1ta^timer^wi^in~^lS~~recorZ 

■ — time  intervals — ■ 


Flush  buffer 
CLOSE  MSFILE 


Write  corrected  MS  record 


etum  to  (0,0 


< 

gesired  MS  record  exist 

yes 

'  _  _ 

c 

Read  desired  MS  record  | 

A. 2  Special  Programs 


Threee  routines  were  developed  to  aid  in  checking  the  files.  The 
first  is  currently  useful  in  checking  data  problems  while  the  latter  two 
routines  are  useful  in  checking  the  contents  of  the  files. 

DUMPMG  prints  records  from  the  archive  tapes.  The  data  is  printed 
on  octal  just  as  it  is  received  by  the  CDC  6600  and  must  be  interpreted 
manually.  This  allows  data  problem  areas  to  be  checked  in  detail. 

PFPRNT  prints  the  contents  of  the  file  AVGEN.  The  number  of  CDC 
words  in  each  record  is  given.  The  label  and  access  data  records  are 
listed.  For  each  data  record,  the  times  and  corresponding  data  are  listed 
in  tabular  form.  The  time  is  given  in  the  concatinated  form  as  well  as 
milliseconds  since  Jan.  0,  1970. 

DMPMS  prints  a  summary  of  the  data  in  the  file  MSFILE.  The  label 
record  is  listed.  For  each  record  containing  data,  the  range  of  times 
covered  by  that  record  are  listed.  The  data  itself  is  not  checked  so  that 
gaps  within  the  time  intervals  are  not  noted. 

Further  development  might  include  the  incorporation  of  these  routines 
into  MGOVL  and  expanding  their  capabilities.  It  would  be  helpful  if  AVGEN 
and  MSFILE  could  be  scanned  in  detail,  within  particular  times  or  in  total, 
and  all  data  gaps  larger  than  a  particular  value  be  reported  for  all 
stations  or  for  selected  ones. 


/OVERLAY/ 


Defines  KEYWORD  -  Overlay  Relation 


NOVL  10 

Number  of  overlavs 

HOVL  16 

Max.  number  of  overlavs 

NDOVL  5 

First  dimension  of  KOVL 

KOVL  (5,16) 

Keyword  overlay  array,  1st  index 

/ SOVRLY/ 

Sequential  Overlay  Link 

KSTK  0 

Index  for  KSTK  for  next  overlay  after 
decriinent 

MSTK  8 

Max.  number  overlay  links 

NSTK  (8) 

Overlay  link  stack 

/INITOV/ 

Initialization  Switch 

INITV  0 

First  card  is  read  as  title  if  INITV 
is  0 

not 

I/O  Common 

IN 

Input  file 

I  OUT 

Output  file 

/I  OX/ 

I/O  Extra  Files 

I  OLD 

Old  AVGEN  to  be  copied  from 

I  AG 

General  average  data  file 

MSF 

Mass  storage  file 

/MGDATC/ 

Data  Input  File 

MGDATA  6LMGDATA 

Magnetometer  Archive  Data  Tape 

/KEY/ 

Keyword  common 

NZ 

Number  of  items 

MZ  12 

Max.  number  of  items 

KEY 

Keyword 

TZ  (11) 

Other  items  from  card 

KKEY 

Keyword  key 

K7  (11) 

Keys  indicating  formats  used  for 
interpreting  IZ 

/ EQUIV/ 

Equivalence  Common 

NEQ  3 

Number  of  defined  parameters 

MEQ  8 

Max.  number  of  parameters 

IEQ (8)  0 

Values  obtained  from  card 

NKQ(2,8) 

Symbol  def ini t ion/expec t ed  key 

TABLE  2.2.0C 

COMMONS  IN  BASE  OVERLAY 
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/CARD/ 

MCARD  80 
KARD(80) 

/TITLEC/ 

DATE (2) 
TITLE(8) 


STA  numerical  (real)  Station  string 
DEL  numerical  (real)  Output  data  step 
ASEC  numerical  (integer)  Number  seconds 

average 

Keyword  card  common 

Number  of  columns  used 
Card  characters 

Run  Title  Common 

Dace  and  time  when  title  was  read 
Title  from  title  card 


/PRTST/ 

PRTEST 

/STLBLC/ 

STLBL(8) 


Print  Switch 
Printing  to  be  done? 
Storage  Label  Common 
Label  for  data  set 


/TIMEC/ 


INTC  0 
INTDAY  * 
INTIME  * 
INTMS  * 
LSTC 

LSTDAY  * 
LSTIME  * 
LSTMS  * 


Title  Common  (*means  initialized  by 
MSTIME) 

Initial  time  concatinated 
Days  since  0  Jan.  1970 
Milliseconds  in  day 
Milliseconds  since  0  Jan.  1970 
Last  time  concatinated 
Days  since  0  Jan.  1970 
Milliseconds  in  day 
Milliseconds  since  0  Jan.  1970 


/STATC/ 


ND  3 

NQT 

KSTA 

JSTA 

NSTA 

MSTA  10 
LBLOST (10)  D 


Station  Common  (*means  initialized  by 
STATNS  D  means  initialzed  by 
data  statement) 

Number  of  data/quantities/station/time 
point 

Number  of  quantities  in  time  point 

Station  sequence  number 

Data  station  number 

Number  of  stations 

Max.  number  of  stations 

Label  external  index 
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IDOST(IO) 

D 

Data  index  (external  index) 

LBLDST(IO) 

Label  (data  index) 

IODST(IO) 

1 

External  index  (data  index) 

ISOST(IO) 

* 

Storage  index  (external  index) 

ISDST(IO) 

* 

Storage  index  (data  index) 

IDLSST(IO) 

* 

Label  (storage  index) 

IOSST(IO) 

External  index  (storage  index) 

IDSST(IO) 

Data  index  (storage  index) 

/LIOBUF/ 

Length  of  Arrays  Stored  in  AVGEN 

LIOAVG 

69 

Length  of  /AVGC/ 

LIOST 

96 

Length  of  /STATC/ 

/PLOTC/ 

Plot  Parameters 

IPLOT 

-1 

Type  of  plot 

MMM 

0 

Plot  Parameter 

MAX 

1 

Plot  Parameter 

NP 

720 

Number  of  joints  plotted 

START 

Start  time 

PLOTIM 

24 

Plot  time  (hr) 

MSPLOT 

Plot  time  (ms) 

DEL 

2. 

Point  spacing  (min) 

MS  DEL 

Point  spacing  (ms) 

DELTIC 

2. 

Time  between  tics  (hr) 

/PLTLBL/ 

Label  for  Plot 

PLOTID(3) 

Label  set  by  data  statement 

/ST/ 

CRT  Plot  Initialization  Parameter 

SF 

1.5 

Plotting  scale  factor 

XP 

.5 

X  -  plotting  frame  value 

YP 

.5 

Y  -  plotting  frame  value 

/INTERC/ 

Common  for  Interactive  Use  Pa-ameter 

KINTER 

0 

If  non-zero,  adjust  for  iuttia  ' i vc 
operation . 
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/DATAK/ 


Buffer  to  Handle  Data  Quantities 


NDATA 

30 

Length  of  data 

DATA( 30) 

Data  buffer 

/OOPSC/ 

Buffer  for  problem  data 

MS  OK 

1500 

Acceptable  millisecond  change  interval 

NDOOPS 

Storage  needed  for  one  data  set  plot 

information 

IOOPS 

Next  error  index  to  be  used 

LOOPS 

96 

Error  buffer  length 

KOOPS(96) 

Error  buffer 

/BUFX/ 

Quantities  Defining  Buffer  Status 

MS  BASE 

-1 

First  time  in  buffer  (ms) ,  used  to 

calculate  index 

MSPT 

1000 

Milliseconds/time  point 

MS  SET 

Milliseconds/data  set 

MPSET 

10 

Number  of  time  points/data  set 

NQS 

3 

Number  of  data  quantities/station  point 

NQST 

Number  of  quantities/output  set 

NQSET 

Number  of  quantities/input  set 

/FFG/ 

Error  Flags 

FFLAG 

65536. 

Fill  flag  when  there  is  no  data 

XFLAG 

32768 

Data  greater  than  this  value  not  used 

AVGFLAG 

65536. 

Value  used  if  there  are  no  points  in 

average 

/AVGC/ 

Quantities  Defining  Averaging 

NMSAVG  =MSTP*NTAVG 

Milliseconds  averaged  over 

NTAVG  10 

Time  points  averaged  over 

MSDAVG  =MSTP*NTDAVG 

Millisecond  shift  between  averages 

NTDAVG  1 0 

Time  point  shift  between  averages 

LAVBUF  =NQST*NTAVG 

Storage  needed  for  one  set  of  averages 

IDAVBF  =NQST*NTDAVG 

Index  range  for  one  set  of  averages 

MN 

Min  index  for  quantities  being  averaged 

MX 

Max  index  for  quantities  being  averaged 

NS2 

60 

Composite  dimension  of  S  and  XN 

S  (  30) 

Sums 

XN ( 30) 

Number  of  items  in  each  sum 

['ABLE  2.2.2C  ADDITIONAL  COMMONS  IN  OVERLAY  (2.2) 
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/mgdimf/ 

See  MGREAD  for  exp]  .mat i on .  NR EC 
set  to  0  for  proper  restart 

// 

Large  buffers 

IOB(512) 

Final  output  buffer 

IBUF(7680) 

Buffer  to  hold  data  before  averaging 

/ TOBUFK/ 

Quantities  Defining  Buffers 

JIOB 

I0B  storage  index 

LIOB 

512 

I0B  length 

1SX 

IBUF  storage  index 

LBUF 

7680 

IBUF  length 

LIMTSX 

IBUF  storage  limit  causing  averaging 
when  crossed  by  initial  item  of  a  set 

NISX 

=0 

Number  of  times  TSX  crossed  L1MISX 

>11 SX  • 

5 

Wax.  number  of  times  LIW1SX  may  be 
crossed  before  mandatory  flushing  of 
buff  er . 

/BMBFIL/ 

Data  to  be  Printed  at  Bombout 

LMS 

CONSTANTS 

Last  time  processed  in  milliseconds 

IDT 

34 

Data  type  wanted  from  MGREAD 

MSMDIF 

200000000 

Max.  forward  time  jump  processed  in 

milliseconds . 
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// 

HDZGAM 

XT  I  ME 
DYGAM 
BSLV 

/FFC/ 

FFLAG 

XFLAG 

AVGFLAG 


/AVGC/ 


NMSAVG 

NTAVG  . 

MSDAVG 

NTDAVG 

LAVBUF 

IDAVBF 

MN 

MX 

NS2 

S(30) 

XN(30) 

/PLOTMS/ 

INPC 

INPDAY 

TNPTIM 

INPMSM 

LSPC 

LSFDAY 

LSPTIM 

LSPMSM 

JPTC 

JPTDAY 

JPTIME 

JPTMS 


TABLE  2.2 


(720,3) 

(720) 

(3,7) 

(3,7) 


65536. 

32768. 

65536. 


10 

10 
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Plotting  Arrays 

Components  of  magnetic  field  at 
different  times 
Corresponding  times 
Half-range  sizes 
Mid-range  values 

Error  Flags  (Same  as  in  overlay  (2,0) 

Fill  flag  when  there  is  no  data 
Data  greater  than  this  value  is  not  used 
Value  used  if  there  are  no  points  in 
average 

Quantities  Defining  Averaging  (As  in 
overlay  (2,0) 

Milliseconds  averaged  over 
Time  points  averaged  over 
Millisecond  shift  between  averages 
Time  point  shift  between  averages 
Storage  needed  for  one  set  of  averages 
Index  range  for  one  set  of  averages 
Min.  index  for  quantities  being  averaged 
Max.  index  for  quantities  being  averaged 
Composite  dimension  of  S  and  XN 
Sums 

Number  of  items  in  each  sum 
Plot  Times  Common 

Initial  time,  concatinated ,  from  data 
set  of  card 
,  day 

,  ms  within  day 
,  ms  since  0  Jan.,  1970 
Last  time,  concatinated,  from  MSXTRA  of 
card 
,  day 

,  ms  within  day 
,  ms  since  0  Jan.,  1970 
Ongoing  plot  time,  concatinated 
,  day 

,  ms  within  day 
,  ms  since  0  Jan.,  1970 


.  3C 


ADDITIONAL  COMMONS  IN  OVERLAY  (3,0) 
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/PLTERR/ 


Error  Check  Common 


PLTXIN 

65536. 

Defines  errors 

PLTXOU 

99999. 

Value  used  for  errors 

/INBUFC/ 

Input  Buffer  Common 

MIOB 

512 

Buffer  size 

IOB (512) 

Buffer 

/ AGLBLC/ 

Avgen  Label  Common 

LAGLBL 

20 

Length 

IAGLBL(20) 

Label 

/RED3WC/ 

Lindage  Needed  for  Restart 

MVIRG 
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/FFC/ 

FFLAG 

XFLAG 

AVGFLAG 


/AVGC/ 


NMSAVG 

NTAVG 

MSDAVG 

NTDAVG 

LAVBUF 

IDAVBF 

MN 

MX 

NS  2 

S  ( 30) 

XN ( 30) 

/ RAMSX/ 

1NTCM 
1NTMSM 
LSTCM 
LSTMSM 
NSDM  (3) 
NSTAM  (7) 

I  MOST 

NOSETM  (21) 
NT  REG 
MSDMS 
MSREC 

ERE  i'  MS  (3  780) 
E RAMSX 

ERE CMS  (3781) 
FT  LI  .MS 


TABLE  2 . 2 . AG 


65536. 

32768. 

65536. 


10 

10 


oO 


180 

120000 


44 

65536. 


Error  Flags  (same  as  In  overlay  (2,0) 

Fill  Flag  when  there  Is  no  data 
Data  greater  than  this  value  not  used 
Value  used  if  there  are  no  points  in 
average 

Quantities  defining  averaging  (as  in 
overlay  (2,0) 

Milliseconds  averaged  over 
Time  points  averaged  over 
Millisecond  shift  between  averages 
Time  point  shift  between  averages 
Storage  needed  for  one  set  of  averages 
Index  range  for  one  set  of  averages 
Min  index  for  quantities  being  averaged 
Max  index  for  quantities  being  averaged 
Composite  dimension  of  S  and  XN 
Sums 

Number  of  items  in  each  sum 

Parameters  Defining  Random  Mass  Storage 

Initial  time  concatenated 
Initial  time  milliseconds 
Final  time  concatinated 
Final  time  milliseconds 
Data  items  per  stations 
Number  of  stations 

Conversion  of  standard  index  to  mass 
storage  index 

Number  of  quantities  per  time  set 
Number  of  time  sets  per  record 
Milliseconds  per  time  set 
Milliseconds  per  record 
Length  of  mass  storage  record  file 
Length  of  saved  common 
Length  of  mass  storage  record 
Fill  value 


ADDITIONAL  COMMONS  IN  OVERLAY  (4  ,0) 
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SEARCH 

NOW 

KREAL 

PRTITL 

IEQUIV 

NUMCVT 

KEYVTRD 

HOLLER 


TABLE  2. 


1 


Searches  for  proper  keyword. 

Obtains  date  and  time  in  characters. 

Creates  real  and  integer  numbers  from 
input  character  string. 

Prints  title,  time,  and  page  numbers. 

Scans  input  for  equal  signs  and 
substitutes . 

Converts  between  real  and  integer  if 
necessary . 

Reads  and  interprets  free  field  input. 

Creates  Hollerith  "values"  from  input 
character  string. 


2.1R 


ROUTINES  IN  OVERLAY  (1,0) 


73 


! 

i 


c 


FILDTA 


MSTOTAL 
MGREAD 
I SEARCH 
ADDFLUX 

MOVE IT 

MGUNPK 

ISLOTER 

ITISUM 

CHECK 

TIMSTOP 

STODTA 

AVGOUT 

IOBOUT 

BMBOUT 

IFCLC8 

NOW 


Extracts  data  from  data  base  and  stores 
it. 

Converts  time  to  milliseconds. 

Extracts  desired  data. 

Positions  data  base. 

Combines  flux  data  which  is  stored 
separately  in  base. 

Adjusts  data  in  double  buffers. 

Unpacks  packed  data. 

Finds  next  data  set. 

Extracts  time  from  data. 

Checks  data  error  flags. 

Adjusts  time  parameters  to  handle  time 
skips . 

Stores  desired  data  in  buffer. 

Averages  data,  stores  average. 

Outputs  averaged  data  and  adjusts  buffer. 
Is  called  if  system  terminates  program. 
Finds  I/O  file  list. 

Obtains  date  and  time  in  characters. 


TABLE  2.2.2R 


ROUT”  j  :N  OVERLAY  (,2,0) 


MGOVL 

STATNS 

MSTIME 

ICTIME 

MOVE 

INTGRS 

IORDR2 


TABLE  2 


1 


Main  routine,  calls  overlays. 

Initializes  station  index  tables. 

Converts  concatinated  time  to 
milliseconds . 

Converts  milliseconds  to  concatinated 
t  ime . 

Efficiently  moves  blocks  of  data. 

Separates  number  into  digits. 

Orders  a  list  and  adjusts  associated 
list . 


2. OR 


ROUTINES  IN  BASE  OVERLAY 


5 


// 


Buffers 


I OB  (512) 

Avgen  buffers 

MSX  (1536) 

Msfile 

index 

MSB  (4096) 

Msfile 

buffer 

/MSBUFK/ 

Buffer 

access  parameters 

L10B 

Current 

length  of  avgen  buffer 

LMSX 

Current 

length  of  msfile 

index 

LMSB 

Current 

length  of  msfile 

buffer 

JIOB 

Current 

index  of  avgen  buffer 

JMSX 

Current 

index  of  amsfile 

index 

JMSB 

Current 

index  of  amsfile 

buffer 

/MSLBLK/ 

Initial 

label  for  msfile 

MSLBL  (10) 

Initial 

label  for  msfile 

SPECGRM 

IRED3W 

BASE 


MSPLOT 

AAXIS 

DAMOYR 

LINER 


Controls  plotting  of  data  from  AVGEN. 

Extracts  data  to  be  plotted. 

Scans  data  and  sets  normalization 
parameters . 

Plots  three  magnetic  quantities. 
Writes  axis  numbers 

Labels  plot  from  concatinated  time. 
Plots  data. 


TABLE  2.2.3R 


ROUTINES  IN  OVERLAY  (3,0) 


FILCVT 

INITMS 

INITAG 

NOW 


TABLE  2. 


Converts  data  between  files. 
Initializes  mass  storage  file,  MSFILE. 
Initializes  AVGEN. 

Obtains  date  and  time  in  characters. 


2.4R 


ROUTINES  IN  OVERLAY  (4,0) 
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3. 


Program  Use 

This  program,  MGOVL,  is  designed  to  be  run  either  interactively 
or  in  batch  mode.  The  program  is  structured  such  that  the  loader  under 
the  CDC  NOSBE  system  will  create  a  disk  file  with  a  local  name  MGOVL. 
This  disk  file  can  be  made  permanent.  The  program  has  been  restricted 
to  less  than  600000  words  of  memory  in  order  to  be  run  interactively. 

The  running  time  may  be  rather  long,  650  sec,  if  entire  tapes  are  pro¬ 
cessed.  The  usual  tasks  done  interactively  use  much  less  time  although 
it  is  desirable  to  raise  the  default  interactive  time  limit. 

A  description  of  the  input  is  given  in  section  3.1.  The  current 
options  are  listed  in  section  3.2. 
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The  input  described  below  is  a  free  field  input  which  is  used  for 
both  batch  and  interactive  use.  A  card  in  batch  input  is  equivalent  to  a 
line  in  interactive  use. 

The  first  card  read  is  saved  as  a  title  card  and  used  to  label  the 
first  page  of  output.  The  remaining  cards  are  either  keyword  cards  or  cards 
which  follow  the  keyword  card  requiring  them.  Each  keyword  causes  values  to 
be  stored  or  an  action  to  be  performed.  The  main  inter-dependence  is  the 
use  by  one  keyword  of  values  set  by  an  earlier  keyword.  The  program  will 
stop  upon  reading  an  END  card  or  an  end  of  record. 

The  keyword  cards  are  free  field  cards.  The  fields  are  separated  by 
the  break  characters  and  Leading  and  trailing  blanks  are 

ignored.  Scanning  ends  when  a  termination  character,  either  ")"  or 
plus  20  blanks,  or  the  end  of  the  card  is  found. 

A  field  with  just  digits  and  possibly  a  preceeding  sign  is  treated 
as  an  integer.  A  real  number  occurs  in  a  field  with  a  point  and/or  an  E 
preceeding  an  integer  field.  All  other  fields  are  alphanumeric  and  their 
first  ten  characters  are  used  starting  with  the  first  non-blank.  Blanks 
will  not  be  removed  from  within  the  field  but  will  be  added  to  the  end  if 
necessary.  After  a  termination  character,  the  rest  of  the  card  can  be  used 
for  comments.  A  null  field,  no  non-blank  character  between  two  breaks  or 
beyond  the  termination,  usually  results  in  defaults  being  used.  The 
fields  bracketing  an  "="  may  form  an  equivalence  pair.  If  the  first  field 
can  be  identified,  the  second  field  is  used  and  both  fields  are  removed 
from  the  card  image.  If  the  first  field  cannot  be  identified,  acts 

as  a  standard  break  character.  If  the  value  after  the  "="  is  numeric,  the 
program  will  convert  it  t~  the  proper  form.  This  is  not  true  in  general. 

The  following  are  the  current  equivalence  keys: 


STA  =  concatinated  list  of  standard  indices  for  the 
stations  to  be  processed.  The  default  is  all 
stations . 

DEL  =  minutes  between  output  data  points  or  averages. 

The  default  is  1/6  minutes  or  10  seconds. 

ASEC  =  seconds  to  be  averaged  over.  The  default  is 
10  seconds. 

If  a  keyword  cannot  be  identified,  an  error  message  is  printed  and 
the  next  card  is  read.  If  the  first  column  of  the  card  contains  a  "C"  and 
no  keyword  can  be  identified,  it  is  assumed  to  be  a  comment  card  labelling 
the  input. 
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3.2  Keyword  IV.  f i nit  ions 

The  present  keywords  and  their  fields  are  defined  on  the  following 
pages.  If  the  card  is  terminated  before  all  fields  are  defined,  those  not 
defined  will  be  treated  as  null  fields.  The  keyword  is  given  in  upper 
case;  the  following  fields  in  lower  case. 

The  concatinated  time  used  for  many  keywords  is  an  integer  in  the 
form:  YYDDDHHMMSS  with  YY,  years  since  1970;  DDD,  Julian  day  within  year; 
HHMMSS,  hours,  minutes,  and  seconds  within  day. 
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TITLE 

NOPRINT 

PRINT 

END 

STATION 

TIME 

STORLBL 

FILEDATA 

ADDATA 

PLOTLBL 

INTER 

N INTER 

SF 

PORG  (x,; 


Reads  next  card  as  a  new  title. 

Stops  printing  of  input  cards. 

Starts  printing  of  input  cards.  This  is  the  default. 

Stops  execution.  An  end-of-record  or  end-of-file  will 
also  halt  the  program. 

(Station  index,  station  index,  station  index ,... station  index). 

Initializes  station  list  based  on  standard  indices;  see 
Table  1.1.  This  function  is  also  done  by  equivalence  key 
"STA' ;  see  section  3.1. 

(Initial  concatinated  time,  final  concatinated  time)  Sets  time 
parameter  for  processing). 

Reads  next  card  as  a  new  title  for  AVGEN  file.  The  default 
title  is  "MAGNETOMETER  NETWORK  DATA.  FIRST  RECORD  CONTAINS 
PACKING  AND  AVERAGING  DATA." 

(Initial  time-concatinated ,  final  tlme-concat inated) 

Initializes  AVGEN  and  outputs  data  found  on  MGDATA  between 
the  times  requested  and  from  the  stations  requested  and 
averaged  as  specified.  See  the  equivalence  keys  in  section 
3.1. 


(Initial  time-concatinated,  final  time-concatinated) 

Positions  AVGEN  to  allow  averaged  data  to  be  transferred 
to  it  from  MGDATA  starting  at  the  initial  time.  The  status 
and  averaging  parameters  are  determined  by  those  previously 
used  for  AVGEN. 

Reads  the  first  thrity  columns  of  the  next  card  to  be  the 
plot  identifier.  The  default  is  "FOUGERE  X3827  MAGNETOGRAM". 

Causes  the  program  to  stop  before  and  after  each  plot  for 
interactive  control. 

Causes  the  program  to  continue  through  the  plots.  This  is 
the  default. 

(Scale  factor)  Sets  the  plotting  scale  factor.  If  null,  it  is 
reset  to  default  1.5. 

Sets  the  plotting  origin.  The  default  is  0.5, 0.5. 

A  null  value  in  either  field  leaves  the  former  value 
unchanged . 
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I 


A. 1  Debug  Disk  File  Structure  and  Use 

The  following  sets  of  control  cards  set  up  a  disk  file  containing 
a  loadable  binary  file  and  an  UPDATE  file.  In  each  set  the  cards  referring 
to  disk  files  must  be  changed  from  application  to  application  and  the  cards 
defining  and/or  saving  the  I/O  file  sets  must  be  added. 

I.  Debug  Initialization 

This  set  of  cards  reads  in  the  Fortran  or  compass  cards 
using  UPDATE.  An  UPDATE  file  is  created  and  the  decks 
are  compiled.  If  either  causes  an  error  the  run  is 
aborted.  If  compilation  is  successful,  a  binary  file  is 
created.  These  files  are  combined  and  put  in  a  permanent 
disk  file  with  the  binary  as  the  first  file.  An  execution 
is  then  attempted. 


-JOB  CARD- 

COMMENT.  DEBUG  (AFGL/NOS)  -  INITIALIZATION 
UPDATE (N,W,) 

FTN  (I=COMPILE,  SL,  T,  R  =3) 

REWIND  (LGO) 

REQUEST  (NL,*PF) 

COPYBF(LGO.NL) 

COPYBF( NEWPL.NL) 

CATALOG (NL,MGOVLX,ID=  ) 

ATTACH  (PLTLIB .TEKLIB) 

ATTACH  (HTWLIB , ID=  ) 

LI BRARY (HTWLIB , PLTLIB) 

MAP (FULL) 

LGO. 

7/8/9 

UPDATE  creation  deck 

7/8/9 

Program  test  data 

7/8/9 

6/7/8/9 
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II.  Debug  Testing 


This  set  of  cards  compiles  changes  to  UPDATE  file  and  sub¬ 
stitutes  them  using  COPYL.  This  is  for  testing  only  and  does 
not  modify  the  disk  file. 


-JOB  CARD- 

COMMENT,  DEBUG  (AFGL/NOS)  -  TESTING 
ATTACH  (OL.MGOVLX,  ID=  ) 

COPYBF (OL , BIN) 

UPDATE  (P=OL,R=C,Q) 

FTN  (I=COMPILE,  SL,  T,  0PT=0,R=3) 
REWIND (BIN, LGO) 

COPY (BIN, LGO, B) 

ATTACH (PLTLIB , TEKLIB ) 

ATTACHE  HTWLIB , ID=  ) 

LIBRARY (HTWLIB .PLTLIB) 

B. 

7/8/9 

UPDATE  Modification  deck 
7/8/9 

Program  test  data. 


III.  Debug  Update 

This  set  of  cards  modifies  the  UPDATE  FILE  and  compiles  a  new 
binary  file.  After  cataloging  the  new  disk  file  a  run  is 
attempted. 


-JOB  CARD- 

COMMENT,  DEBUG  (AFGL/NOS  - 
ATTACH  (OL.MGOVLX  ID= 
SKIPF(0L,1,17) 

UPDATE  (N,P=OL,R=CN,W,F) 

FTN ( I=COMP ILE , SL , T , R=  3 ) 
REWIND (LGO) 

REQUEST (NL,*PF) 

COPYBF (LGO, NL) 

COPYBF (NEWPL.NL) 

CATALOG (NL .MGOVLX , ID= 

ATTACH (PLTLIB .TEKLIB) 

ATTACH (HTWLIB, ID= 

L I BRARY ( HTWLIB , PLTL I B ) 

MAP (FULL) 

LGO. 

7/8/9 

UPDATE  modification  deck. 


UPDATE 

) 


) 

) 


(continued  on  next  page) 
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7/8/9 

Program  test  deck 
7/8/9 
5/7 / 8/9 

IV.  Production  Run 

This  set  of  cards  directly  loads  the  program. 

-JOB  CARD- 

ATTACH ( OL , MGOVLX , ID=  ) 

ATTACH (PLTLIB , TEKLIB ) 

ATTACH(HTWLIB,ID=  ) 

LIBRARY(HTWLIB, PLTLIB) 

OL. 

7/8/9 

Data  Deck 
7/8/9 
6/1/8/9 

V.  Debug  Audit 

These  control  cardw  ill  list  the  binary  decks  in  the  binary 
file  and  the  decks  in  the  UPDATE  FILE  and  will  compile  all 
decks  to  produce  a  listing.  The  list  of  binary  decks  in¬ 
cludes  the  date  of  compilation,  entry-point,  and  subroutines 
of  each  routine. 


-JOB  CARD- 

COMMENT,  DEBUG  (AFGL/NOS)  -  AUDIT 
ATTACH(OL, MGOVLX, ID=  ) 

DISPOSE (OUTPUT , *PR=C) 

SKIPF (OL ,1,17) 

UPDATE (P=OL , R=C , F , L=F) 
FIN(I=COMPILE,SL,T,R=3) 

REWIND (OL) 

ATTACH ( PLTL I B , TEKLIB ) 

ATTACH (HTWLIB, ID  =  ) 

L I B RARY ( HTWLI B , PLTL I B ) 

MAP  (FULL) 

LOAD(OL) 

NOGO. 

6/ 7 / 8/9 
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SOME  PROPERTIES  OF  SOLUTIONS  OF  THE  PHOTOELECTRON  FLUX  EQUATIONS 
The  solution  of 


-ay"(x)  -  by'(xO  +  cy(x)  -  dy(x+A)  =  S(x) 


on  (0,oo)  which 


1. )  is  bounded  on  [0,co), 

2. )  is  continuous  on  f0,oo),  and 

3. )  has  continuous  first  derivative  y’(x)  on  (0,cd)  with  the 

property  that  lim  y'(x)  exists. 

x— »  0 

Let  S  =  u  +  iv  be  a  complex  variable  and  let 
K  (S)  =  S2  -  q  +  PeAS 


with 


4ac~Hb 

4a2 


and  P  = 


2a 


There  can  be  at  most  one  real  zero  of  K  (S)  which,  if  it  occurs, 

must  occur  at  So  =  — - - f~~  +  ^ '  This  zero  may  be  of  order  1  or 

A2 

order  2.  All  other  zeros  of  K  (S)  are  of  order  1  and  occur  in  conjugate 
pairs.  Let  +  iv^  for  k  =  1,  ....  n  be  all  the  zeros  of  K  (S)  to  the 
left  of  the  line  Re(s)  =  b^2a  and  \  —  ivk  ior  ^  =  n+l,...be  t  lie  zeros  of 
K(S)  to  the  right  of  Re(s)  =  b 1 2a' 

(There -are  a  countably  infinite  number  of  these). 


CASE  I. 


-x/ 

S(x)  =  Sq  e  xq  +  S^r (x)  where  r(x)  is  a  rectangular  pulse  of 


V 


width  2r  and  height  1/2T  centered  at  x 


Let 


1 


H(x)  = 


(1/ 


-x/xo  +  (*.+»)  «  (  A  -  £ 


X  2 


-A/, 


o  —q  +  pe  X 


n  111  h 

A  i/2  +  q)  X  +  S  (i^-  ~) 


k=l 


Cos  (v  +  i  ^ ) 


1A.  For  X>X  +  T, 


m  /  , 

V(x)  -  Hte)  +  E  Bk«  !“k  -  -j;  <X-(X1-D  )  Cos  + 

-  K=1  k  k 


+  ck  e  (uk  -  IT  >  (x-(xi+r)  >  Cos  (vk  X  +  V  H  k; 


IB. 


For  Xx  -  T'fXtXj  +  Y , 


Y(x)  =  H (x)  +  X  B  e  k  2a 
k=l  k 


(u,  -  ^  )  (x  -  (x -D 


1  *  '  Cos  (  kX  +  W  ) 


1C.  For  0<X<X1  -  r  » 


E  C  e  (uk  2a  *  Cx-(xx+r)  >  Cos  (v  +  Q  ) 
k=n+l  k  k  k 


Y(x)  =  I!(x)  +  E 

h=n+l 


B  e(uk  ~  2a  )  (X"(X1  +  T)  )  Cos  (v  X  +  W,  ) 
k  k  k 


+  C  e  <Uk  ”  2a  )  <X-(Xi+r)  >  cos  (  v  X  +  6  ) 
k  k  k 


If  S  is  not  a  zero  of  K  (S)  then  A=B=0 
o  - 

If  S  is  a  zero  of  K  (S)  of  order  one  then  A=0  but  B  is  arbitrary, 
o 

If  S  is  a  zero  of  K  (s)  of  order  two  then  A  and  B  are  both  arbitrary, 
o 
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2k.  For  X>XL, 

Y(x)  =  H(x)  +'  X  BKe  (V  2a  )  (X~V  Cos  (vkX+ 

2B.  For  O^X  <X, 

Y(x)  =  H)x_  +  2  B  8  (uk  Ta  J  (X_X1)  Cos  (v.  +  ) 

k=m+l  K  k 


All  the  remarks  about  the  constants,  location  of  zeros  of  K(s),  etc 
that  followed  Case  I  also  can  be  applied  to  this  case. 


-  DETAILS  - 


The  solution  of 


-ay"  (x)  -  by'  (xO  +  cy  (x)  -  d  y  (x+A)  ■  S(x) 


on  [O,00) 

which 

1.) 

is  bounded  on  [0,co), 

2.) 

is  continuous  on  [0,ao),  and 

3.) 

has  continuous  first  derivative  y'(x)  on 

lim  y'  (x)  exists 

(0,oo)  such  that 

x — >0 

y"(x)  +  ;  y  (x)  -  ;  y(x)  +  ~  y  (x-A)  - 

9  cl  ot 

-  a  S(x) 

Change  the  dependent  variable  from  y  to  Z  via  the  relation 

..  b 

y  =  Z  (x)  e  a  X‘ 

This  yields  the  equation 

(//)  Z  "  (xO  -  q  Z  (x)  +  p  Z  (x+A)  -  S*(x) 
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where 


q 


0,  p 


b 

2a  A 


0, 


* 


S  (x) 


b 

2a 


x  S(x) 


If  y(x)  has  properties  1,  2,  and  3  listed  above  then  Z  (xy)  is 
(e  ^  x)  and  has  properties  2  and  3.  Conversely,  if  Z  (x)  is 

O'  (e  ^  x)  and  has  properties  2  and  3  then  Y(x)  has  properties  1,  2,  and  3. 

If  Z,  (x)  is  a  solution  of  the  homogenous  equation 
(*)  Z  (x)  -  aZ(x)  +  pZ(x+A)  =  0 

which  is  0  (e  *)  and  has  properties  2  and  3  and  Z^(x)  is  a  solution  of 
#  with  the  same  properties,  then  Z^  (x)  +  Z^  (x)  is  also  a  solution  of  // 
with  these  same  properties. 


If  we  assume  that  Z(x)  is  a  solution  of  the  homogenous  equation  * 
which  is  0(e  ^  x)  and  has  properties  2  and  3,  then  the  Laplace  transform 

/co 

Z  (x)  e  sx  dx 

0 


converges  on  Re  (S)  ^  and  transforming  all  of  *  yields 

A 

Z^o)  S  +  Z  (0)  e  SA  f  Z  (t)  e  "St 
Z  (s)  -  +  KTsT  q 


dt 


where  K(s) 


q  +  pe 


SA 


The  zeros  of  K  (s)  . 

Setting  S  =  u+iv  in  -  q  +  pe^  =  0,  separating  the  real  and 
imaginary  parts  and  setting  them  equal  to  zero,  and  doing  a  little  more 
algebra  one  finds  that  the  zeros  of  K(s)  are  in  the  intersection  of  the 
curves 
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The  curves  of  u  and  V  are  plotted  for  a  sufficiently  large  value  of 
p.  The  .curves  are  symmetric  w.r.t.  the  u  axis  so  the  lower  halves  have  not 


been  plotted. 


K'  (s)  “  2s  +  pAe^ 

2  SA 

K"  (sO  -  2  +  pA  e 

K  (s)  -  K'  (s)  -  0  implies  that  S  -  ^  + 

Therefore  there  is  only  one  possible  zero  of  order  two  of  K.  (s)  at 

S  =  i 
o  A2 

K"  (s)  =  0  implies  that  V  =  ^  so  there  can  be  no  zeros  of 

order  three  or  more. 

All  zeros  of  K  (s)  are  1st  order  except  perhaps  a  zero  of  order  two 

at  S  . 
o 

All  zeros  (except  perhaps  the  zero  at  So)  occur  in  conjugate  pairs. 
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There  are  only  a  finite  number  of  zeros  to  the  left  of  Re(s)  =  z~ 


There  are  c  ountably  infinite  number  of  zeros  to  the  right  of  Re(s)  =  j" 


If  p  is  sufficiently  small  there  will  be  no  zeros  of  K  (s)  to  the 


left  of  Re(s)  *  r-  . 

za 


It  is  extremely  unlikely  that  a  zero  of  K  (s)  lies  on  the  line 
Re(s)  =  ^  (probability  zero),  but  it  is  possible. 


If  there  are  no  zeros  of  K  (s)  on  Re(s)  =  —  one  can  use  the 

b 

Mellin  inversion  integral  taken  along  the  contour  Re(s)  =  to  find  the 
solution  of  the  homogeneous  equation  *  . 

(.)  ,  /£  +  i0°  (Z(o)  S  +  z'(o)  e  SXjds 

(s)  2  i  ^  - — 


L 


it  _  |G0 
2a 


2ni 


b  +  j  oo 

— 

/2l 

SA 

e 

£  Z  (t)  eSt  dt 

/ 

K(s)  , 

Jo 

/  b 

- 

2a  -  i 


co 


St  , 

e  ds . 


For  X?0  one  may  complete  contours  for  the  first  integral  using 
semi-circles  in  the  half  plane  Re(s)  £  ^  and  use  the  Cauchy  residue 
theorem  to  find  that 

b 


2a 


+  i 


oo 


1 

2ni 


(Z  (o)  S+z’_(oJi  efj.  .  (Ax+B)e(  I  ~  Ta  '  aI  +  q) 


S  2  q  +  peSA 


2a 


-  i 


OO 


™  S,  x.  -  S.X 
+  E  aRe  k  +  ake  k 

K=  1 


=  (Ax+B)e 


( i  -  —  -  /rr 

v  A  2a  \/A2 


q) 


3u  u,  ^  Cos  (v^  X  + 


+  E  A.  e  k 

92  K=1  ^ 


k  i 


We  now  turn  to  the  evaluation  of 


2a  +  1  °° 


S(x=(A-t)) 


ds  converges  uniformly  for  te  [0,S],  hence 


2.  _  j  00 

2a 


«  z  (t)/2"* + 1  >  ds 


.  i 

2a 


K(s) 


/, 


2^  +  1  00 


—  .  i  00 
2a 


S 

K(s) 


f. 


Z(t)  e  ~St  dt 


Sx  . 
e  dx. 


(and  A-  t  >  0  for  te  [0,A]  ) 


so  as  long  as  X  >  0  one  may  again  complete  contours  using  semi¬ 
circles  in  the  half  plane  Re(s)  <  and  the  Cauchy  residue  theorem  to 

find  that 


1 

2ni 


+  too 


7 


—  -  1  00 
2a  1 


S(x  +  (At)) 


K(s) 


-  (Ax+B)  e 


M 


SA  (x+(A-t)) 


+  E  b  e  Sk 

K*=l  k 


sk(X+bke) 
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2rti  /  dtZ<c> 


+  i  o° 


ii  -  i  OO 

2a 


S (x+(A-t) )  „  , 

- K(7j -  ds  =  (A*  x  +  B’)  e  o 


,  “  ,  S.  X  b*  eSk  X 

•  2j  0  .  6  K  T  K 

K=1  k 


s  x  m  u 

(A'x+B')  e  o  +  2  Be 

K=1 


Cos  (  VkA  +  WR) 


Since  it  is  easily  verified  that  e  k  is  a  solution  of  the 

homogenous  equation  *  if  and  only  if  S  is  a  zero  or  K(s),  and  in  case  S 

S  X  k  ° 

is  a  zero  of  order  2  than  Xe  o  is  a  solution  of  *,  it  follows  that 

the  constrants  A,  B,  A^,  A',  B',  and  B^  in  the  above  evaluation  of  the  two 

integrals  are  all  arbitrary.  We  may  combine  the  results  of  both  integrations 

to  find  that  the  general  solution  of  *  which  is  &  (e  ■  X)  and  satisfied 

properties  2  and  3  is, 

S  X  A  _  \  X  Cos  (  v,  +  <t> ) 

(Ax  +  B)  e  o  +  E  k  k 

K=1  ^ 

» 

S,  =  u.  +  i  v,  and  S,  are  the  zeros  to  the  left  of  Re(s)  =  r-  except  for 
k  k  k  k  2a 

possibly  a  zero  at  S  .  If  S  is  not  a  zero  of  K  (s)  then  A=B=0.  If  S  is 

oo  o 

a  zero  of  order  1,  then  A=0. 


If  one  now  assumes  that  Z(x)  is  a  solution  of  the  non-homogeneous 
equation  (X  which  is  O'  (e  X)  and  satisfies  2  and  3,  then  taking  the 
Laplace  transform  of  //  and  inverting  leads  to; 


Even  though  Z(x)  is  now  a  solution  to  the  non-homogeneous  equation  if 
rather  than  the  homogeneous  equation  *,  the  evaluation  of  the  first  two 
integrals  proceeds  in  exactly  the  same  way  as  in  the  homogeneous  case  and 
gives  only  solutions  to  the  homogeneous  equation  *.  Therefore,  the  most 
general  solution  of 

Z"  (x)  -  q  Z  (x)  +  p  Z  (x+A)  =  S*  (x) 


of  d  (e  rr  x)  which  satisfied  properties  2  and  3,  provided  it  exists,  is 
Zb 


given  by 


m 


Z(x)  =  (Ax+B)e  SoX  +  £  A^e  ^  Cos  (VkX  +  ^k) 

K=1 


^  +  i  oo  *  ,  s 

1  /2a  s+  Cs) 

2i  /  — i ir/TT- 


K(s) 


e  x  ds . 


CD 
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One  can  now  check  Z(x)  ar  i  find  that  it  does  satisfy  it  and  is  of  ^  (e  X) , 
and  satisfied  properties  2  and  3. 

The  results  fcr  Case  I  and  Case  II  are  obtained  by  evaluating  the  last 
integral  in  the  expression  for  Z(x)  and  multiplying  Z(x)  by  e  X).  Keep 
in  mind  that  for  certain  values  of  X  contours  must  be  completed  in  the  half 
plane  Re(s)  <  b/,2a  while  for  other  values  of  X  the  contours  must  be  com¬ 
pleted  in  the  half  plane  (Re(s)  2a .  This  last  procedure  leads  to  the 

sum  of  an  infinite  number  of  terms  because  there  are  an  infinite  number  of 
roots  of  K(s)  in  the  right  half  plane. 


Solutions  to  special  cases  of 

co 


(1)  C(x)Z(x)  =  S(x)  +  f  K(r,x) 

_/x+A 


Z(r)  dr. 


In  the  paper  "Degradation  and  High  Energy  Scattering  of  Radiation"  V. 
Fano  solves  certain  equations  of  the  type 

co 

U(E)  n  (E)  =  f  (E+r ,r)  n  (E+r)  dr  +  S(E) 

•/) 

by  observing  that  when  S(E)  is  the  unit  impulse  function  S(E-Eo)  the  equation 
can,  under  certain  restrictions  on  the  kernel,  be  converted  to  an  equation 
in  which  the  integral  is  a  convolution  integral.  This  makes  the  resulting 
equation  amenable  to  transform  techniques.  If  S(x)  is  the  unit  impulse 
function  equation  (1)  can  also  be  converted  to  an  integral  equation  in 
which,  again  under  certain  conditions,  the  integral  can  be  written  as  a 
convolution  integral.  The  resulting  equation  still  contains  a  term  with  an 
advanced  argument,  but  in  many  cases  this  equation  can  also  be  solved  using 
Laplace  transforms. 


First,  assume  that  c(x)  >■  0  on  [0.°°)  an<i  let  y(x) 
k(r,x)  =  K(r,x)/C(r).  Then 


(1A)  Y(x) 


S(x) 


x+A 


k(r,x)  y  (r)  dr 


c(x)Z(x)  and 
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Otic 


Also  assume  that  S(x)  *  S(x-X  )  and  that  X  >2A.  Define  a 

o  o 

u(x)  in  three  parts  as  follows: 


A.  -  on  OtX<X  -  2A,  u(x)  is  a  solution  of  the  equation 


u(x) 

J  x+A 


k(r,x)  u(r)  dr  +  k(X  ,x) 

o 


B.  -  on  X  -  2  AcXCX  -A,  u(x)  «  k(X  ,x) 

n  *  o  o 


C.  -  on  X>X  -A,  u(x)  =  0 
o 


Now  let  f(x)  *  u(x)  +  S(x-Xq).  Then 


oo 


oo 


oo 


f  k(r,x)  f  (r)  dr  *  /*  k(r,x)  u  (r)  +  /^  kCr.x)^  (r 
-^f+A  *^c+A  x+A 


X  -A 


A.  0<X<X  -  2A 

-  o 


X+A 


k(r,x)  u  (r)  dr  +  k(X  ,x)  =  u(x) 

o 


B.  X  -  2A <X<X  -A,  k  (X  ,x)  -  u(x) 
0—0  o 


C.  X>X  -  A,  0  *  u(x) 
—  o 


Thus, 


J  k(r ,x)  f  (r)  dr  +  cf  (x-x0)  *  u(x)  +  cf  (x-Xq) 


X+ 


function 


-X  )  dr 
o 


f(x) 
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Therefore,  what  one  really  needs  to  solve  (1A)  in  this  case  (S(x)  • 
S(x-Xq))  is  the  solution  of  the  integral  equation 


(2) 


U(x) 


X  -A 


k(r,x)  u(r)  dr  +  k(X  ,x) 

o 


on  the  interval  0  <  X  <X  -2A.  However,  note  that  r  varies  up  to  X  -A 
—  o  o 

in  the  above  integral  so  that  the  values  of  u  (r)  up  to  X^  -  enter  into 
equation  (2) . 

To  cast  equation  (2)  in  a  different  form,  first  change  the  variable 

of  integration  ly  letting  s=X  -A  -r.  Then  r=X  -A  -s  and  0  <  s  C  X  -2A-X. 

o  o  —  —  o 

Equation  (2)  becomes 


(3)  u(x)  = 


X  -X  -2A 
o 


k(XQ  -A  -s,x)  u  (Xq  -A  -s)  ds  +  K(Xo,x) 


observe  that  X  -A  -s  varies  from  X  -A  to  x+A. 
o  o 
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Now  change  the  independent  variable  by  letting  t  =  Xq  -X  -2A.  Then 
X  **  X  -2A  -t  and,  since  0  O  (  X  -2A  we  find  that  0  C  t  ^  X  -2A. 
Equation  (3)  now  becomes 


(4) 


u  (X  -2A  -t) 
o 


k  (X  -A  -s,  X  -2At)  u  (X  -A  -s)  ds 

O  O  O' 

+  kCX  ,X  -2A  -t) 
o  o 


Note  that  even  though  X  -2A  -t  varies  from  0  to  X  -2A,  X  -A- 

o  o  o 

varies  up  to  X  -A. 

o 
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Now  change  the  unknown  function  from  u  to  v  by  letting 
u(t)  =  v  (Xq  -A  -t).  Equation  (4)  becomes 

t 


(5)  v  (t+s) 


k  (X  -A  -s,  X  -2  -t)  (s)  ds+k  (X  ,X  -2A  -t) 

o  o  o  o 


That  is,  if  one  can  find  a  solution  of  (5)  for  0  <  t  <X  -2A  (which 


means  also  that  A  £  t  +  A  ^  X  -A)  then  by  letting  u(t)  -  v(X  -A  -t)  we 


find  that  u  (f)  satisfies  equation  (4), 


Initially,  one  is  interested  in  functions  v  (t)  defined  on 


-A]  which  satisfy  (5.)  for  t  in  [0,X^  -2A] .  One  could  try  to  solve 


(5.)  by  generalizing  standard  integral  techniques  (Neumann  series).  This 
possibility  will  be  investigated  later.  Instead,  we  assume  that 


k  (Xq  -A  -  ,  Xq  -2A  -t)  =  g  (t-s) 


over  the  region  t>0  and  0  <  s  C  t,  and  that  the  Laplace  transform  of  g(x) 
exists.  An  example  of  such  a  kerhel  will  be  given  later.  Under  the  assump 
tion,  equation  (5)  becomes 


(6)  v  ( t  -*-A)  = 


■/ 


g  (t-s)  v  (s)  ds  +  g  (t+A) 


We  look  for  solutions  to  (6)  over  f 0 , 00 ) .  These  will  also  be  solu¬ 
tions  over  [0,  -A] .  (It  is  also  possible  to  look  for  solutions  to  (6) 

on  [0,  X^  -A]  using  finite  Fourier  transforms  and  trying  to  work  out  a  con¬ 
volution  integral  theory  for  such  transforms.  This  possibility  will  also 
be  investigated  later) . 

One  should  not  forget  that  on  [X  -2  ,  X  -A]  we  must  have  u(t)  = 

o  o 

K(X  ,t).  Therefore,  for  0  <■  t  <  A, 


Now  take  the  Laplace  transform  of  equation  (6.) 


V  (t  +A)  e  pt  dt 


G(p)  V(p)  + 


g  (t  +A)  e  pt  dt 


Let  W  *  t  +  A.  Then  we  have  (Here  G(p)  and  V  (p)  are  the  transforms  of  g(x) 
and  (x) ,  respectively. 


v(w)e  'P(W‘S)  dw 


QD 


G(p)  V(p)  +  f  g(w)e 

-/a 


~P (w~a) 


dw 


G(p)  V (p)  _  epA 


g(w)e 


g(w)e  pW  dw. 


e?  V(p)  =  G(p)  V(p)  +  e  PA  G(p) 


*  eP*  G(p) 


It  can  be  shown  that  the  zeros  of  epA  -G(p)  lie  in  a  half  plane  Rep£a 
for  a  sufficiently  large.  If  one  can  also  find  an  a  such  that  the  inte¬ 
gral  in  ati°o 


G-(Pl  e  pt  dp 
epA  -  G(p) 
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converges  uniformly  for  t>o  and  r — is  of  order  0  |  p  |  *)  for 

epA  -G(p)  I  I 

some  k>0,  then  (8.)  defines  a  function  on  [0,°°)  which  satisfies  both  equa¬ 
tion  (6.)  and  boundary  (7.). 

To  show  the  function  defined  in  *8)  satisfies  equation  (6),  we  proceed 
as  follows: 


J  8  (t"s)  (s)  ds  "  2ni  / 

S  n  A  — 1 


.  a+ioo  PA.  . 

/  <p  is^ 

i-ioo  e  -G(p) 


g(t-s)  e  pS  ds 


a+ioo  aAC 

i  /  «  pA  G<p)  .  p*  ( 

2,11  /  epa  -®(p)  / 


g  (d)  e  1X1  da 


ra-iao 


a+: 

=  4/ 

-'a-ioo 


a+ioo 


i  / 

!ixi  y 

“'a— 1 


a+ioo 


e  v  G(p)  e  pt 
ep  -G(p) 


e  pi  G<p)2  '  Pt 
ep4  -G(p> 


lG(P>-/ 


8  (a)  e 


-pd  dal 


dp  -  I  g  (a)  do. 


a+ioo 


pA  x  P(ta)  , 
e  K  G(p)  dp 

-ioo  pA  N 
e  -G(p) 


If  — - (p) -  is  0  (  /p/^)  for  some  R>  0  then  the  contour  in  the 

e  pA  -G(p) 
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last  integral  in  the  above  equation  can  be  completed  to  the  right.  Since 
e  pA  -G(p)  has  no  zeros  in  this  half  plane 


i  r 

2tU  / 


a+ioo 


a-ioo 


£  G(p)  .  P(t-a) 
ePA  -  G(p) 


dp  =  0  if  t  -  a<0 


Thus, 

t 


L 


g  (t-s)  v (s)  ds  -  2^ 


e  PA  G(p) 


( 


G(p) 


e  PA  -  G(p) 


-) 


e  pC  dp 


a+ioo 


1 

2ni 


a-ioo 


e  PA  G(p) 


pA 


e  pA  -  G(p) 


-1 


e  pt  dp 


a+ioo 


1. 

’2ni 


r -Loo 


e  PA  G(p)  e  P  (t+A) 
e  PA  -  G(p) 


a+i 


dP  -  of 


1 

2ui 


G  (p)  e  P(t+A)  dp 


'a-i 


=  Y  (t  +  A)  -  g  (t  +  A) 


/ 


Hence,  v(t  +  A)  =  /  g  (t-s)  v  (s)  ds  +  g  (t+s) 

■'o 


To  show  that  v  (t)  *  g(t)  for  0  £t  <_s,  observe  that 
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v(t) 


/a+ico 

a-ioo 


f.PA-P.(P). 

e  PA  -G(p) 


e  pt  dp 


Jra+io°  ?  \ 

(G  <p)  +  SPA  —  ■  )* Pi dp 
.-loo  V  " 


_1 

2ni 


a-ioo 

■/ a-ioo 


G  (p)  e  pt  dp  +  ± 


a+ioo 

^  f 

■/a-ioo 


£j£L_  e  Pt  dp 
PA  -  G(p) 


But 


a+ioo 

(t)  +  ±  f  9..(SL)..JL.£(J  -A).  dp 

-4-ioo  l‘e  “P  G<P> 


-  a+ioo 

7 

■'a-ioo 


G  (p)2  ep  (t~  ) 


-  e  "P  G(p) 


dp  -  0  if  0  t  <L  because  in 


this  case  the  contour  may  be  completed  in  the  right  half  plane  Re  a  >  o  and 
1-e  P^  G(p)  has  no  zeros  there  fend  G(p)  is  analytic  there). 
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g(t) 


a+i 

1  r 

2  Til  1 

-'a-i 


Pt 


P+1 


dp  provided  Re  a  >  -1 


(9)  V  (t)  -  ±, 


i  ra+1 

-Xi 


p(t+A) 


Cp+1)  epA-l 


dp 


The  zeros  of  (p+i)  e  -1  are  located  on  the  intersection  of  the 


curves 


and 


X  «  -(1+y  cot  y  s) 


2-2  2 

y  «  e  -  (x+ir 


and  all  lie  to  the  left  of  the  imaginary  axis.  They  also  occur  in 
conjugate  pairs. 


Thus,  in  the  integral  in  (9.)  we  take  a  so  that  Z  (o  where  Z 

o 

is  the  real  part  of  the  first  zeros  to  the  left  of  the  imaginary  axis. 


H(p)  = 


(p+l)  e  pA  -1 


is  of  0  (  /p/  ) 


as  p-)co  in  the  half  plane  Re  p  >  a  and  H(p)  is  analytic  there  so  v(t) 
defined  in  (9)  is  a  solution  of  (6)  in  this  special  case  provided  the  in¬ 
tegral  in  (9)  converges  uniformly.  Since  H(p)  is  0  (  /p/  h  as  p-*oo  in 
the  left  half  plane  also,  we  can  complete  the  contour  of  the  integral  in  (9) 
in  the  left  half  plane  also  and  find  that  (assuming  all  poles  are  of  order 


A(pt  +  1)  e 


i=l  A  -pi 


CO 

where  |pi|  are  t^ie  roots  of  (p+l)  e  -1  =  0. 

Summing  by  pairing  the  complex  conjugates,  we  find  that 

N  2e  XiC  [(A-XJ  cos  Y  v  „ 

v(t)  =  Z  L  1  1  t4yl  sin  Yi  c 


2  2 
(A-X.)  +  Y± 


where  p^  =  X^  +  iY^.  This  last  series  is  uniformly  convergent  for  t  >  0, 
so  (9)  is  also  uniformly  convergent  to  t>0. 


In  this  special  case,  the  solution  to  (1  A)  is 


ANGULAR  BEHAVIOR  OF  PHOTOPRODUCTION 


The  upper  atmosphere  has  many  effects  of  interest,  amont  them  are  re¬ 
fraction  and  attenuation  of  radio  waves,  absorption  of  ultraviolet  radiation, 
and  the  slowing  of  low  satellites.  The  behavior  of  the  upper  atmosphere  is 
controlled  by  the  sun  through  many  links,  primarily  the  ultraviolet  radiation 
whose  effect  varies  diurnaly  through  the  solar  zenith  angle.  Other  links 
which  may  be  significant  but  which  are  harder  to  model  are  through  solar 
particles  which  are  funnelled  by  the  earth's  magnetic  field  and  may  be 
approximately  diurnal,  as  the  solar  wind,  or  sporadic,  as  particle  bursts 
from  solar  flares.  A  detailed  model  is  needed  to  understand  the  relative 
importance  of  different  causes  and  effects  and  to  identify  the  limits  of 
various  approximations  and  empirical  relations. 

To  model  the  upper  atmosphere  in  detail  is  difficult  because  it 
varies  with  time  and  with  position,  distance  above  the  earth.  An  even 
greater  and  more  significant  problem  is  the  fact  that  the  ionosphere  is  a 
plasma  which  requires  the  modification  of  the  Boltzmann  equation,  the  stand¬ 
ard  basic  equation  of  kinetic  theory,  through  the  use  of  Fokker-Planck 
methods.  Due  to  the  nonlinearty  of  the  resulting  integro-dif ferential  equa¬ 
tion  no  solutions  are  known  except  for  very  simple  cases.  Dr.  Jasper se  has 
approached  this  modelling  problem  by  first  modelling  only  the  electron  dis¬ 
tribution  function  using  as  given  a  neutral  model  atmosphere,  the  solar 
radiation,  and  the  exospheric  temperature.  The  energy  distribution  of  the 
atoms,  molecules,  and  ions  were  assumed  to  be  Maxwellian.  The  resulting  non¬ 
linear  integro-dif ferential  equation  had  terms  for  photoionization,  recom¬ 
bination,  electron-neutral  elastic  collisions,  excitation  collisions,  de¬ 
excitation  collisions,  impact  ionization  collisions,  and  collisions  between 
charged  particles.  A  stable  iterative  calculational  method  was  found  to 
solve  the  equation  numerically.  The  program  has  been  used  to  model  the  atmos¬ 
phere  measured  by  a  rocket  flight.  Further  checking  is  desirable.  This 
work  should  be  continued  to  identify  the  significant  interactions  and  parame¬ 
ters  in  the  present  model. 
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With  the  understanding  gained,  further  work  would  consider  effects 
not  included  in  the  present  model.  One  such  is  diffusion  but  its  inclusion 
would  require  more  detailed  calculational  modelling  of  the  vertical  structure 
of  the  atmosphere.  The  inclusion  of  diffusion  and  of  modelling  the  density 
of  the  atmospheric  constituents  would  allow  one  to  study  and  possibly  predict 
Che  effect  of  solar  activity  on  radio  communication  and  the  slowing  of 
satellites. 


LO 


Let  the  initial  state,  I,  with  spin  projection  M  be  represented  by 
1 1  M  >  .  Let  |  F'  M* ,  pu> be  the  final  state  of  the  ion  in  State  F'  with 
spin  projection  M'  and  an  electron  with  momentum  p  and  spin  projection  u. 
The  interaction  with  a  photon  of  momentum  k  and  polarization  e  is 


where 


and 


S 


t  £  ik.r 
A  =  e  e 


A  ^  i  V  r 

H  =  A  x  A  =ukne 


/.  A 

n  =  k  ?  e 


The  cross  section  is 


da 
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da  , 
M  Mu 


21+1 


E 

MM 


'M? ,p  u/O/IM^ 


2 


Using  the  relations  from  Brink  and  Satchler’s  book  "Angular  Momentum" 
on  page  151 : 


T.  -A  1  1*1** 

6  (P-P)  =  4Tt-n  2  6  (P-P’)  S  <£  (p)  (p) 

hti  p  ,  m  m 

lm 

e  Lk*r  =  S  (-i)k  (2K+1)  3/2  Jk  (kr)  ^Ck  (k)  Ck  (y) 
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where 


\e  V 

C*  (k)  CK  (r) 


2  Ck  (k)  Ck,  (r)  (km  ki'  /O  0  > 
•  id  m  * 

mm 


The  final  state  and  the  interaction  operator  can  be  expanded: 

<F*  M\  pu>  =  /  F'  M'  p/m  |  u  ;  FM  >  A  F/^(p)  (?)^  * 

m 


<y=  2  |  (°ek  (  {k)  k_1  *)k  )  °  +  ^Mk  (ook  kJ  ° 


^  2  P  \  P\  n  ^ 

where  (p)  =  . '  — -  C  m  (p) 

m  (2 i)  ! 

is  a  symmetric  traceless  sensor  formed  from  the  vector  p  and  0^  and 

k  k  2 

0^  are  the  properly  normalized  elective  and  magnetic  2—  pole  operators. 

Usint  the  reduced  matrix  elements  the  sums  over  the  projections  can 
be  carried  out: 

,  „  1/4  1  k  A  Ak'-1  'O  k'  )  ° 

do  =  2  (factors^)  (  (p)  ^ (k)  U  )  (p)  ^  (k)  Tt J  J 

A  /\  a  a 

where  n  and  n.'  are  e  and/or  H  .  The  summations  can  be  rearranged  via  6-j 

and  9-j  symbols : 

do  =  2  (factors.,)  (( (?)  ^  (p)  ^(k)k-1njk  (((k)k_1n*) k')  J 

=  2  (factors3)  ((p)  ^  (((k)k-1(k)k’-1)  £)*)  )  • 

Summing  over  the  photon  polarization  directions 

2(n  n')  a  (k)^ 

A 

since  k  is  the  only  remaining  direction. 
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Thus 


da  »  E  (factors^)  |(p)^  ^(k)^  (k)™)^)° 

=  E  (factors^)  (  (p)  ^  (k)  ^ J  ° 

=  E  (factors)  P-  (cos  0  ) 


where 

A  a 
cos  &  =  p  .  k 

If  the  highest  multiple  considered  is  2*cinax,  since  /=  k  ^x)  k' , 

0  0kmax 

X  max  =  2 


Parity  is  conserved  and  the  states  have  well  determined  parity. 


The  non-zero  operators  must  all  have  a  definite  parity.  Therefore,  all 
A  k-l  A  k 

terms  (  (k)  n  )  of  interest  must  have  the  same  parity  and 


.  A  A  lr 

<(  (k)kitx)k 
reduced  to  (k) 


a  k'-l*  k'  X 

(  (  k  )  Tt’)  )  must  have  positive  parity.  This  is 
which  has  parity  (-)-'.  Thus  JV.  must  be  even.  A  detailed 


reduction  leads  to  the  same  result. 
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ATMOSPHERIC  TURBULENCE  MODELLING  IMPACTING 
STRATOSPHERIC  ENVIRONMENT  PROGRAM 


Efforts  conducted  under  this  research  impact  the  stratosphere  chemi¬ 
cal — Transport  Model  Development.  The  objectives  and  relevance  are  stated 
as  follows: 

OBJECTIVES :  Develop  methods  of  describing  the  constituents  and  thermal 

distributions  of  the  normal  stratopshere  and  calculate  how  these  distribu¬ 
tions  are  affected  by  Air  Force  flight  operations  and  atmospheric  transport 
coefficients;  characterize  stratospheric  disturbances  that  can  be  produced 
from  mountain  waves  and  underlying  thunderstorm  activity.  Develop  special¬ 
ized  numerical  techniques  to  incorporate  radiation  transfer  and  couple 
photochemical  processes  with  the  heat  balance  of  the  stratosphere. 

RELEVANCE:  This  task  is  the  key  to  the  stratospheric  environment  project 

wherein  the  results  of  the  various  measurements  are  used  to  establish  a 
realistic  model  of  the  stratosphere.  The  model,  in  turn,  is  used  to  predict 
if  environmental  changes  will  occur  as  a  result  of  AF  flight  operation.  This 
task  is  essential  to  the  preparation  of  encironmental  impact  statements  of 
the  Air  Force. 

The  major  effort  of  this  work  here  on  turbulence  mode] ling  is  to  im¬ 
prove  stratospheric  models  for  the  prediction  of  Air  Force  systems  operations 
in  the  stratosphere.  In  particular,  the  probability  of  occurrence  of 
turbulence  for  a  given  location  and  season  as  a  function  of  altitude  is  a 
required  input  to  the  tropospheric/stratospheric  model  used  to  predict 
effects  of  Air  Force  Systems  emission  pollutants  in  the  stratosphere.  Thus 
we  are  investigating  the  seasonal  variability  of  the  characteristic  layers 
of  turbulence  ia  the  troposphere  and  stratosphere  and  determining  the  season¬ 
al  variations  in  the  natural  levels  of  NO  trace  constituents  which  are  also 

x 

engine  combustion  by-products.  This  will  result  in  refinements  in  the 
troposphere/stratosphere  1 -D  and  3-D  chemical/radiative  transport  models 


which  are  used  to  predict  environmental  effects,  such  as  ozone  depletion, 
caused  by  Air  Force  Systems  emissions  in  the  troposphere  and  stratosphere. 


Meteorological  data  (temperature,  winds,  pressure  and  relative 
humidity)  in  the  height  range  0-60  km  on  magnetic  tape  is  currently  being 
restructured  to  provide  efficient,  easy  access  for  processing  with  minimum 
computer  time.  A  statistical  approach  is  being  used  to  determine  the 
seasonal  probability  of  occurrence  of  turbulence  as  a  function  of  latitude, 
longitude  and  altitude  through  the  use  of  the  Richardson  criteria  for  the 
presence  or  absence  of  turbulence  in  the  atmosphere.  N0^  trace  constituents 
measured  by  LKD  will  be  compared  with  data  from  other  experimental  tech¬ 
niques  to  determine  the  extremes  of  seasonal  variability  in  these  constitu¬ 
ents.  The  turbulence  and  NO^  seasonal  results  will  be  used  to  parameterize 
the  important  role  of  transport  processes  and  provide  fine  adjustments  in 
the  trapospheric/stratospheric  chemistry  of  the  1-D  and  3-D  models. 

Introduction: 

The  stratospheric  seasonal  variations  task  under  the  Stratospheric 
Environment  Project  is  currently  studying  the  general  aspects  of  turbulence 
in  the  troposphere  and  lower  stratosphere.  The  data  being  used  in  this 
study  has  been  provided  by  the  Environmental  Data  Service  National  Climatic 
Center  in  Asheville,  North  Carolina.  The  "rawinsonde  data"  consists  of 
winds,  temperature  and  pressure  as  a  function  of  altitude  obtained  from  all 
available  stations  (currently  144  stations)  for  the  period  1948-1976.  At 
the  present  time  only  1970-1976  data  have  been  reformated  for  processing  on 
the  CDC  6600  computer  system.  This  data  is  being  used  in  calculating  the 
Richardson  number,  a  stability  criteria  for  determining  the  presence  or 
absence  of  atmospheric  turbulence,  and  is  given  by  the  following  relation- 


where:  g  is  the  acceleration  of  gravity 


3T/_  is  the  vertical  temperature  gradient 

dz 

T  is  the  dry  adiabatic  lapse  rate  (  10°  K/jq,,) 

andf  fl]2  "  [_  (3V/x/  3  z  )  2  +  (3Vy/  sz5  2]  18  the  Square 

of  the  vertical  shear  of  the  horizontal  winds. 


The  rawinsonde  system  is  not  ideally  suited  to  take  measure¬ 
ments  to  provide  thermodynamic  data  at  the  required  sampling  intervals 
for  this  application.  The  data  being  used  in  this  study  was  intended 
to  provide  northern  hemisphere  pressure  maps  for  a  number  of  milibar 
pressure  levels.  In  order  to  determine  the  Richardson  number  at  equal 
height  intervals,  the  temperature  and  wind  component  data  pointB  are 
fit  using  a  Hermite  interpolation  algorithm.  The  temperature  function 
along  with  the  derivatives  of  the  winds  and  temperature  with  respect 
to  altitude  are  interpolated  at  1/10  km  intervals  and  used  to  calcu¬ 
late  the  Richardson  number.  The  probability  of  occurrence  of  turbulence 
is  obtained  by  defining  a  critical  Richardson  number  (R^  -=  1  to  1/4 
according  to  current  literature)  and  assuming  turbulence  present  when 
the  Richardson  number  is  equal  to  or  less  than  this  value.  The 
probability  is  then  taken  as  the  ratio  of  the  total  number  of  occurrences 
to  the  total  number  of  measurements  at  a  1.0  km  level  over  a  season. 

These  seasonal  results  are  then  averaged  over  one  km  and  plotted  as  in 
the  figure  provided.  Two  techniques  are  used  to  establish  the  seasonal 
tropopause  height.  The  seasonally  averaged  temperature  gradient  value 
is  tested  for  a  sustained  change  in  sign  which  indicates  the  change  in 
lapse  rate  from  the  troposphere  to  the  stratosphere.  Frequently,  at 
high  latitudes,  the  temperature  remains  fairly  constant  above  the 
troposphere  and  through  the  stratosphere.  Since  a  negative  lapse  rate 
may  not  develop,  a  second  test  is  used  which  established  the  tropopause 
height  as  the  maximum  rate  of  change  in  the  temperature  gradient. 

Initial  results  from  stations  at  low,  middle  and  high  latitudes 
reveal  interesting  seasonal  and  latitudinal  variation  in  the  occurrence 
of  turbulence  as  well  as  particular  altitude  levels  which  consistently 
reveal  relatively  high  occurrences  of  turbulence.  It  is  also  interesting 
to  note  that  turbulence  rapidly  decreases  at  heights  which  generally 
agree  with  those  levels  reported  in  models  as  the  tropopause  level  and 
follows  this  pattern  as  the  tropopause  level  decreases  from  its  highest 
level  in  equatorial  regions  to  its  lowest  levels  at  high  latitudes. 


•i 
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This  particular  result  is  substantiated  in  theory  since  the  stratosphere 
region  above  the  tropopause  level  is  a  dry  stable  region  having  a  nega¬ 
tive  temperature  lapse  rate,  and  thus  the  degree  of  turbulent  mixing 
and  transport  is  greatly  reduced  from  that  occurring  in  the  troposphere. 

The  resolution  of  the  temperature  and  wind  data  (  100m  to  a 
few  kilometers)  varies  with  the  rise  rate  of  the  rawinsonde  balloon  and 
the  criteria  of  selecting  data  points  between  mandatory  levels  (those 
standard  tnilibar  height  levels  recorded  at  all  stations) .  The  fit  ap¬ 
plied  to  the  data  to  provide  Richardson  numbers  number  determination 
by  possibly  introducing  effects  not  present  in  the  thermodynamic  structure 
of  the  atmosphere.  The  possibility  of  obtaining  a  limited  amount  of  data 
at  suitable  sampling  intervals  to  obtain  some  knowledge  of  the  actual 
structure  in  the  data  has  been  investigated.  Personnel  at  the  Chatham,  MA. 
rawinsonde  station  have  cooperated  in  the  initial  investigation  by  demon¬ 
strating  the  present  techniques  in  data  acquisition  during  a  visit  to 
that  site  in  September  of  1978.  They  also  provided  at  that  time  an  analog 
strip  chart  recording  of  temperature  along  with  the  standard  printout  of 
thermodynamic  data.  Data  points  extracted  from  the  temperature  analog  at 
the  smallest  possible  interval  indicate  spectrum  wave  lengths  of  approxi¬ 
mately  a  km  and  longer. 


TURBULENCE  STUDY  -  STATISTICAL  ANALYSIS 


Objectives: 

The  statistical  analysis  plan  to  be  described  below  has  two  major 
objectives: 

1)  To  relate  occurrences  of  turbulence  to  latitude,  longitude  and 
seasonal  variation. 

2)  To  identify  areas  of  high  turbulence. 


I  -  Site  Selection 


At  present  data  are  available  on  144  sites  (locations)  for  the  years 
1971-1976.  From  these  sites  we  have  selected  a  subset  of  14  sites 
for  initial  analysis.  This  subset  has  been  selected  to  include: 

1)  A  wide  range  of  latitude  and  longitude.  The  14  sets  are 
representative  of  sites  ranging  from  70°  to  120°  longitude 
(W)  and  25°  to  50°  latitude  (N) .  This  will  give  the  results 
broad  applicability. 

2)  "Possibly"  interesting  variables  such  as  land  configurations 
(i.e.  coastal,  mountains,  flat  lands).  This  will  allow  us  to 
investigate  variables  that  may  introduce  variability  into  our 
results. 

3)  Stations  which  are  geographically  close.  This  will  allow  us 
to  investigate  the  correlation  between  nereby  sites. 

4)  Stations  that  also  have  data  from  1948  to  1970.  This  will  per¬ 
mit  a  more  extensive  time  series  analysis  (if  deemed  useful). 

The  results  obtained  from  the  14  selected  sites  should  be  repre¬ 
sentative  of  locations  in  the  continental  United  States.  When  the  analysis 
is  completed  on  these  sites  data  from  other  sites  will  be  analyzed  to  act 
as  a  validation  of  the  results  obtained  from  the  14  original  sites. 

II  -  Data  for  Analysis 

Dependent  Variable  -  The  major  dependent  variable  is  the  proportion 
of  days  for  a  given  time  period  (e.g.,  a  season)  with  turbulence. 
Turbulence  is  determined  according  to  whether  the  Richardson  number 
is  less  than  .25.  At  each  site  we  have  for  a  given  baloon  launch 
at  most  one  Richardson  value  for  each  1.0k  interval  from  ground  leve 


up  to  a  height  of  35  km.  A  typical  situation  is  to  have  approxi¬ 
mately  two  balloon  launches  per  day  providing  two  readings  for  up 
to  25  km.  At  this  point  the  balloon  usually  bursts  and  measure¬ 
ments  are  no  longer  available. 

Independent  Variable  -  The  independent  variables  include:  site 
latitude,  site  longitude,  altitude,  average  temperature,  North- 
South  wind  speed,  East-West  wind  speed,  temperature  variation 
(e.g.,  range  of  temperature  over  a  season),  density,  humidity 
year,  and  season  (i.e.,  appropriate  time  period). 


Analysis  of  Individual  Sites  (First  Analysis) 

The  first  analysis  is  a  site-by-site  analysis.  That  is,  each  site 
is  analyzed  individually. 

A  -  Summary  Statistics 

For  a  given  time  period  (season)  summary  statistics  on  all 
variables  will  be  computed.  These  include  means ,  standard 
deviations,  maximums/minimums  and  ranges  for  each  altitude 
bins  (each  1km  is  referred  to  as  a  bin) . 

At  the  present  time  we  are  in  the  stage  of  the  analysis.  In 
particular  we  are  determining  the  number  of  altitude  bins 
that  can  be  analyzed  with  good  statistical  accuracy  and  pre¬ 
cision.  High  altitude  bins  (say  above  25  km)  often  do  not 
have  sufficient  data  for  a  good  analysis.  Based  on  the  results 
we  obtain  from  the  14  selected  sites  we  will  determine  how  many 
bins  should  be  carried  further  in  our  analysis. 


B  -  Altitude  Analysis 

For  each  altitude  bin,  and  time  period  (season)  a  two-way 
analysis  of  variance  will  be  performed  using  Time  Period  by 
year  as  the  variables  for  classification.  The  objectives  will 
be:  (1)  quantify  yearly  variation  and  (2)  quantify  time 

period  variation. 


C  -  Altitude  by  Time  Period  by  Years  Analysis 

A  three  way  analysis  of  variance  will  be  performed.  This 
analysis  will  quantify  variation  due  to  altitude,  variation 
due  to  years,  variation  due  to  time  period  and  interactions 
among  these  three  sources. 


D  -  Regression  Analysis 


A  multiple  regression  analysis  will  be  performed.  The  regression 
model  will  relate  proportion  of  turbulence  to  year,  time  period, 
altitude  and  weather  variables  suc.h  as  temperature,  density, 
humidity,  wind  direction,  and  interaction  of  these  variables. 
Specifically,  the  independent  variables  are:  Year ,  Time  Period, 
(season  or  month).  Altitude,  Average  Temperature,  Temperature 
Variation,  Wind  Direction,  (North-South,  East-West),  Density , 
Humidity  and  interaction  of  these. 

IV  -  Combination  of  Sites  (Second  Analysis) 

The  second  analysis  consists  of  combining  sites  and  analyzing  them 
jointly.  If  deemed  necessary  the  following  analysis  can  be  performed  on 
subsets  of  sites. 

A  -  Regression  Analysis  (altitude  x  latitude  x  longitude) 

A  regression  relating  proportion  of  turbulence  to  geographic  location 
and  altitude  will  be  performed.  This  regression  will  also  include 
year  and  time  period  (season) .  The  model  is 

(Proportion  turbulence)  =  b^  +  b^  (Altitude)  +  b^  (latitude)  +  b^ 
(Longitude)  +  (Year)  +  b,.  (time  period). 

If  necexxary  the  model  will  be  expanded  to  include  interactions  of 
the  above  variables. 

B  -  Expanded  Regression  Analysis 

The  next  analysis  will  further  expand  on  the  above  regressions. 
Variables  which  quanify  land  configurations  will  be  added.  Also 
weather  variables  as  discussed  in  III  D  will  be  added. 

V  -  Logistic  Regression 

The  objective  of  analysis  is  to  relate  occurrence  of  turbulence  to 
sets  of  variables.  In  particular  we  want  to  relate  variables  to  the 
probability  of  turbulence.  A  useful  model  for  this  is 


P  = 


1  -t 


1 


E  b .  Xi+b 
e  i  o 
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where  P  represents  the  probability  of  turbulence  and  £  b^  X^  + 
represents  a  linear  function  relating  variables  . . Xfl  to  the  probabili¬ 

ty  P.  By  use  of  the  transformation 


we  have 


b  +  E  b .X. 
o  i  i 


and  regression  analysis  is  now  possible  using  the  transformed  dependent 
,  variable.  We  will  perform  these  logistic  regressions  in  addition  to  the 
avove  regressions. 


A  -  Individual  Sites 

A  logistic  regression  will  be  performed  using  variables  listed 
in  III  I). 

B  -  Combination  of  Sites 

Logistic  Regressions  will  be  performed  using  variables  described 
in  IV  A  and  also  IV  B. 

Note:  Analysis  IV  A,  IV  B,  V  B  will  allow  us  to  relate  the  occurrence 

of  turbulence  to  latitude,  longitude,  altitude  and  seasonal  variation. 
Additionally  these  analyses  will  quantify  importance  of  weather 
variables  to  turbulence.  The  functions  developed  in  these  analyses, 
especially  the  results  of  V  B,  will  permit  us  to  identify  the  profiles 
of  variables  that  relate  to  higher  turbulence  -  viz,  we  will  be  able 
to  identify  values  of  XI,..., X  given  in  V  that  produce  high 
probabilities  P  given  by  the  logistic  regresstion. 

Note:  While  the  above  analyses  address  mainly  the  problem  of  explaining 
past  data,  we  will  also  investigate  their  ability  to  predict 
turbulence.  ' 

Note:  As  stated  earlier  the  above  will  first  be  performed  on  the  14 

selected  sites.  When  this  is  completed  other  sites  in  the  continent¬ 
al  United  States  will  be  analyzed  for  validation.  Then  high  latitude 
sites  will  be  selected  for  ana  analysis  similar  to  the  one  described 
above . 
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PROCESSING  OF  PLASMA  FREQUENCY  PROBE  DATA 


1 .  INTRODUCTION 

Purpose  of  the  Experiment 

The  plasma  frequency  probe  is  a  rocketborne  instrument  used  for 
in-situ  measurements  of  electron  density  and  temperature  in  the  E-  and  F- 
regions  of  the  ionosphere.  This  model  measures  the  parallel  resonant 
frequency  of  an  electrically  short,  RF  excited  antenna  immersed  in  the 
ionospheric  plasma.  Measurements  taken  can  be  used  to  find  electron 
density  and  temperature  and  when  correlated  with  altitude,  an  electron 
density  profile  can  be  derived. 

Theory  of  Operation 

The  plasma  frequency  probe  (PFP)  makes  Its  measurements  by  finding 
the  parallel  resonant  frequency  of  an  RF  excited  antenna. 

The  RF  excitation  is  generated  by  mixing  the  output  of  a  107  to  147 
MHz  sweeping  oscillator  with  the  outputs  of  two  crystal  oscillators 
(106.5  MHz  and  106.4  MHz),  resulting  in  two  signals  sweeping  from  1  to  40 
MHz  precisely  100  kHz  apart.  The  sweep  oscillator  Is  a  varactor-tuned 
oscillator  driven  by  a  ramp  control  voltage,  (sweep  voltage).  This  sweep 
voltage  swings  from  -7.5  to  13  V  in  20  ms  and  then  waits  at  13  V.  As  the 
sweep  voltage  increases,  the  frequency  at  the  oscillator  output  decreases. 
The  two  outputs  are  then  filtered  and  amplified  to  drive  the  antenna.  The 
antenna  current  is  measured  and  this  signal  is  mixed  with  the  second  signal 
from  the  RF  oscillators  down  to  a  100  kHz  signal.  At  parallel  resonance, 
the  antenna  current  is  at  a  minimum  causing  the  output  signal  from  the  RF 
head  to  be  at  a  minimum.  This  minimum  is  then  detected  by  a  peak  detector. 
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The  peak  detector  finds  a  minimum  and  starts  the  frequency  measuring 
sequency.  At  a  minimum  a  flip  flop  is  set,  turning  on  a  counter  in  the 
sequencer.  This  counter  then  counts  1  ms  and  resets  the  flip  flop  in  the 
peak  detector  making  a  1  ms  pulse,  the  timing  gate  pulse,  when  a  minimum 
occurs.  This  timing  gate  pulse  stops  the  sweep  generator,  holding  the 
oscillator  frequency  constant,  while  also  enabling  the  frequency  counter 
for  1  ms.  The  output  of  the  counter  is  the  frequency  in  kHz  (in  binary) 
of  the  minimum.  After  the  sequence  is  completed,  the  peak  detector  is 
disabled  until  it  is  reset  by  the  sequencer  at  the  beginning  of  the  next 
frame . 


At  the  end  of  each  frame,  the  plasma  frequency,  then  stored  in 
counters,  is  jammed  into  a  shift  register.  The  data  is  then  fed  to  the 
buffer  amplifier  serially  during  the  following  frame  at  a  rate  of  500  bits 
per  sec.  The  sequencer  generates  the  reset  and  jam  pulses  along  with  the 
500  Hz  clock  used  by  the  shift  register.  Also,  the  sequencer  counts  the 
1  ras  pulse. 

There  are  three  outputs;  the  reset  pulse,  the  peak  detector  output, 
and  the  plasma  frequency.  Each  of  these  is  fed  through  a  buffer  amplifier 
to  protect  the  telemetry  from  over-voltages  and  at  an  output  impedence  of 
Ik.  To  give  it  the  correct  format,  the  plasma  frequency  data  is  buffered 
by  a  variable  gain  amplifier.  The  amplifier  has  a  gain  of  one  for  the 
reference  bits  (the  first  four  bits),  and  a  gain  of  one-half  for  the 
remaining  bits. 

Probl ems  With  the  Data  and  Results  Obtained 

Figures  1  and  2  show  the  raw  electron  density  data  obtained  from 
the  experiment.  Figure  2  is  an  expanded  view  of  the  interesting  portion 
of  the  flight  upleg  and  downleg.  Due  to  instrument  problems,  a  very  high 
data  drop  out  rate  and  system  noise  level  were  present. 
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Standard  data  reduction  techniques  would  be  useless  with  this 
poor  data.  However,  by  applying  modern  statistical  estimation  techniques 
(weighted  minimum  variance  estimation)  the  electron  density  profile  of 
Figure  3  was  obtained.  The  weighting  was  based  on  the  signal  bandwidth 
and  the  length  of  the  data  set  available. 
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